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ABSTRACT 


Reliably modeling noise attenuation through interaction with vibrating bound- 
ary structures is fundamental to the formulation of effective active noise control sys- 
tems. In this paper we investigate, through numerical approximation, uniform ex- 
ponential stability of two systems which model the acoustic/structure interaction of 
an air-filled, rectangular cavity. The first model assumes dissipative boundary condi- 
tions along one side of the boundary, while the second assumes dissipative boundary 
conditions along all four sides of the cavity. We obtain weak variational formulations 
for these models, express each as finite dimensional systems, and use the Galerkin 
technique to transform the distributed parameter systems into systems of ordinary 
differential equations. We analyze the stability of the finite dimensional systems in 
order to gain insight into the stability of the original infinite dimensional systems. 
Essentially, our analysis consists of solving a generalized eigenvalue problem and ob- 
serving where the eigenvalues lie within the complex plane. This stability analysis 
leads us to conclude that one model is better suited for use in the formulation of the 


noise control problem. Accesion For 
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I. INTRODUCTION 


Imagine that you are ten years old and your teenage brother asks you to 
participate in a scientific experiment. You agree. Inside a 55 gallon drum you go. 
Your brother, clever lad that he is, manages to suspend the drum (and you) a few feet 
in the air. So there you are, dangling by a rope three feet off the ground-sealed in a 
metal drum. With your baseball bat in hand, he strikes the drum dead center. After 
the drum stops reverberating, he lowers it to the ground. As he pops the lid off, he 
asks you to describe how the acoustic noise field in the drum behaved after his forceful 
swing. As a composed budding young scientist, you disregard the blood dripping from 
your left earlobe and answer, “Initially, a large noise field was created by the strike of 
the bat. However, after the impulse force was experienced, the intensity of the noise 
field steadily decreased until it eventually was imperceptible.” Your brother thanks 
you dearly and then proposes a heat conduction experiment... 

In this paper, we too examine noise transmission attenuation through a vi- 
brating boundary structure. Our approach, however, is to do so through numerical 
approximation. More specifically, we examine the exponential stability of infinite di- 
mensional second order systems of coupled partial differential equations by numerical 
approximation. This is a topic attracting considerable interest in the fields of engi- 
neering and applied mathematics because the applications are both numerous and 
diverse—reducing noise levels in automobiles, aircraft, and space launch vehicles to 
name afew. The degree to which a control system formulated to effect noise reduc- 
tion in a fluid-filled cavity succeeds depends to great extent on how accurately the 
underlying mathematical model agrees with the observed physical phenomena. As 
our young scientist reported above, the acoustic field created by an external force 
acting on the cavity boundary steadily diminished to zero as time passed (in the 
absence of any sustaining force). This is equivalent to requiring that all solutions to 


the system of equations selected to model the behavior of the acoustic field, as well 


as the vibrating boundary, converge exponentially to zero by a uniform rate of decay. 
Establishing this result for the infinite dimensional system is nontrivial and, at times, 
a very difficult task. 

To help clarify these ideas, consider the following abstract formulation. Ac- 
cording to [Ref. 1], many examples related to acoustics or fluid/structure interactions 


can be modeled abstractly as a first order system of equations as follows: 
y(t) = Ay(t), t>0, y(t)EH, (1.1) 


where H is an appropriately defined Hilbert space. The companion linear control 


system is typically written 
y(t) = Ay(t)+ BA(t), Ah(t)hEeR", (1.2) 


where A is a control] input and B 1s a linear operator from Rk” into H. This system 
possesses uniform exponential stability if there exist M > 0 and 6 > 0 such that for 
all t > 0 and for all (y(0), y:(0)) € H 


lly(t), we(t)Ilee S Me*"||(y(0), y(0) lu [Ref.2, 3], 


where |ly(t), yz(t)||z, denotes the energy of the system at time ¢ and ||(y(0), y:(0)||z 
denotes the energy of the system at t = 0. 

The authors of [Ref. 1] state that “the most common approach for the ap- 
proximation of a control problem involving I.2 is to formulate a sequence of finite 


dimensional control systems of the form 
ys (t) = ANy"(t)+ BYVAL), t>0, y*() EH", (1.3) 


where the dimension of the finite solution space H” increases toward infinity as N 
tends to infinity. In general, equation I.3 is derived from I.2 using space discretization 
techniques such as finite difference, finite elements or spectral methods developed 
for the approximation of the solutions of I.1. A control strategy is then designed 


for the finite dimensional control problem involving I.3. This control is used as an 











approximation to the desired control function for the infinite dimensional control 
problem I.2. One of the most practical conditions - assure the well-posedness of 
the finite dimensional control problem, as well as the convergence of the approximate 
controls to the desired control for the infinite-dimensional system, is that the solutions 
of 1.3 for h = 0 preserve the exponential decay of the solutions of I.1.” Hence 


determining the stability of the finite dimensional system 
yYHSA Ye), 150, PA OeEH (1.4) 


is of fundamental importance. 

Typically this stability analysis is accomplished by examining where the eigen- 
values of A% lie within the complex plane, since theory tells us that system I.4 is 
globally stable if and only if all eigenvalues of A% have negative real parts [Ref. 4]. 
The process is nontrivial because numerical results have indicated that many popular 
approximation schemes fail to maintain a uniform decay rate as the dimension of the 
approximating system I.4 increases [Ref. 1], even in cases where the original system 
was proven to be exponentially stable. 

In this paper, we examine two different infinite dimensional models by numer- 
ical approximation in order to gain insight into their adequacy for control system 
formulations. Model I is believed to be stable, but not uniformly exponentially sta- 
ble [Ref. 5]. Model II is believed to be exponentially stable [Ref. 6]. 

In the following chapters, topics introduced above are addressed in greater 
- detail. In Chapter II, we develop a general and two specific models describing an 
acoustic field inside a two-dimensional fluid-filled cavity surrounded by a perturbable 
boundary. In Chapter II], we illustrate the Galerkin technique chosen for our numer- 
ical approximations, and in Chapter IV, we apply this technique to the models under 
consideration. In Chapter V, we describe how specific approximations were obtained 
and present our results. We conclude with summary comments and propose future 


areas of study in Chapter VI. 











Il. MATHEMATICAL MODELS 





Force 


Figure 1. 2-D Acoustic Chamber for the General Mathematical Model 


In this chapter, we develop a general and two specific models for a system con- 
sisting of the wave equation coupled with the beam equation on part of the boundary. 
Consider the two-dimensional rectangular air filled cavity surrounded by an impene- 
trable boundary shown in Figure 1. A noise source exterior to the cavity produces a 
perturbing force f which induces vibrations in the cavity boundary causing fluctua- 
tions (i.e., undesirable noise) in the acoustic pressure field within the cavity. Equa- 
tions of motion describing boundary vibrations and acoustic pressure fluctuations 
within the cavity, together with appropriate initial value and boundary conditions 
form a system of coupled, second order partial differential equations in time. 

We begin formulation of the general model by assuming that the interior acous- 
tic pressure field satisfies the standard wave equation ¢ = c?Ag@ where ¢ is the 


velocity potential throughout the cavity and c is the uniform speed of sound in the 





fluid. The velocity potential ¢ is a complex-valued function satisfying 
u(t, L, y) = —V o(t, L, y) ’ 


where ¥ denotes the fluid velocity. Initial value conditions (0,2, y) and ¢,(0,z, y), as 
well as boundary conditions along 0Q—either Dirichlet, Neumann or dissipative—are 
specified. 

A Newtonian analysis of forces and bending moments leads to equations de- 
scribing the motion of the elastic walls bounding the cavity. For simplicity we assume 
that the boundary walls are impenetrable and that only one side of the rectangular 
boundary is perturbable (The term beam refers to the perturbable Bedary,): We 
assume an Huler-Bernoull: beam where (i) w(t,x) denotes the transverse displace- 
ment of the beam of length a, (ii) p, and py denote the uniform mass densities of the 
beam and fluid, respectively, (iii) M(t,2) denotes the total internal moment of the 
beam, and (iv) f(t,2) denotes the external forcing term. The beam equation takes 


the general form 
Powe (t, x) Ba Mee (t, Z) aa —p soit, L, w(t, t)) + TG x) ’ (II.1) 


where —p¢,(t, x, w(t, x)) is the acoustic pressure (This is the first coupling term we 
see in the acoustic/structure system.). 

Initial value conditions are specified, and boundary conditions for the beam 
indicate whether the ends are free, partially restrained, or clamped. Additionally, 
assumptions specifying the types of internal moments of the beam are necessary. 


Internal moments typically consist of bending and damping moments, 
M(t, x) = Mbending SH Maacing ) 


where 


Mbending = Strain component = E(xr)I(x)wz,(t, x). 


The stiffness of the beam is given by E(x)I(x), where E(x) is the Young’s modulus 


and I(x) is the cross-sectional area of the beam. Muamping is taken to be either 














Kelvin-Voigt, spatial hysteresis, or time hysteresis damping. While spatial hysteresis 
damping is an appropriate choice for a beam constructed of composite materials and 
time hysteresis damping is an appropriate choice for a beam constructed of “material 
with memory”, Kelvin-Voigt damping assumes a memoryless beam of uniform (linear) 


mass density where the damping stress is proportional to the strain rate. That is, 
Kelvin-Voigt damping: Miamping = strain rate component = cp{xz)I (xr) weet , 


where cp/(z) is the product of the Kelvin-Voigt damping coefficient cp(x) and beam 
cross-sectional area I(x) [Ref. 7]. Here we consider a beam of uniform cross sec- 
tional area and density. Thus we assume Kelvin-Voigt damping and take the coeffi- 
cient functions cp(zr)I(x) and E(z)I(z) to be constant for all z (i-e., ep(x)I(x) —> 
cpl and E(z)I(r) — ET). 

Based upon the above discussion, the generalized mathematical model for the 


acoustic/structure system is: 


Det = c*Ad for (x,y) E eee t>Q, 
powr(t, ©) + Mzz(t,x) = —prdi(t,z, w(t,z))+ f(t,z), 0<a2<a, t>0, 


(11.2) 


with appropriate initial value and boundary conditions specified for the wave and 
beam equations. 

Next we present two specific models and formally show that each model can 
be expressed in both weak and strong formulations—which, as we shall see, are equiv- 
alent given appropriate choices of inner product spaces. This preliminary work will 
formally justify expressing Models I and II as first order systems in time thus facil- 
itating our numerical examinations of solution stability. In Chapter III we examine 
in some detail the popular variational scheme, the Galerkin method, used to obtain 
approximate solutions to coupled systems such as those addressed in this paper. We 
use the Galerkin scheme to transform the infinite dimensional system II.2 into a finite 


dimensional system to facilitate numerical analysis. 





A. MATHEMATICAL MODEL I 





Force 


Figure 2. 2-D Acoustic Chamber for Model I 


The theoretical model considered in this section is shown in Figure 2. Here the 
wave equation is coupled with the beam equation on one side of a 2-D rectangular, 


air filled cavity. We assume Neumann boundary conditions along [’. That is, 
Vo-n = 0, (z,y) ET, t>0, 


where 7 represents the the outer normal with respect to [. The boundary is, in 
effect, a sound-insulated rigid wall which prevents any acoustic energy from escaping, 
or alternatively, the boundary is a perfect reflector of acoustic waves. 

Further, we assume that the perturbable boundary (i.e., beam) can be charac- 
terized as an impenetrable fixed-end Euler-Bernoulli beam with Kelvin- Voigt damping 
and that both ends of the beam are clamped. Given these assumptions, the equations 


of motion describing the vibrations of the perturbable boundary are: 


powi(t, xc) + Miz(t,z) = —prdi(t,z,w(t,z)) + f(t,z), O<x<a, t>0 (I13) 





w(t,0) = w,(t,0) = w(t,a) = we(t,a) = 0, t>0, (11.4) 


where 


w(t,z) = the transverse displacement of the beam 
o(t,z,y) = the fluid velocity potential 
pp = the linear mass density of the beam 
pf = the uniform mass density of the fluid 
M(t,z) = Elwec + cplweet 


f(t,z) = the force due to an exterior noise field. 


For the stability analysis of this model, we consider the open loop problem absent any 
exterior noise field (1.e., f(¢,z) = 0 for ¢ > 0) and do not concern ourselves with any 
noise control aspects. We assume (i) that the beam is impenetrable to the adjoining 


fluid and obtain the second coupling equation (i.e., the continuity of velocity) 
wi(z,t) = Vod(t,z,w(z,t))-n , O< x<a,t>0, (II.5) 


and (ii) that displacements from the beam’s position of rest are small, which is in- 
herent in the Euler-Bernoulli formulation. Because of (ii), we take the transverse 
displacement of the beam as w(t,z) = wW(t,r) + 6 where © = 0. Under these 


assumptions, equation IJ.3 and II.5 become 
prwr(t, r) a5 Mzz(t, x) = —psOr(t, x, w(t, x) + 6) (11.6) 


w(t,z) = Vod(t,z,w(t,z) + d)-an. (11.7) 


By using two term Taylor Series expansions of ¢; and V¢ with respect to y about the 


point z, equations II.6 and II.7 become 


prwre(t, v) - M(t, x) sk —pslOx(t, x, 0) + Pty (t, x, 0)w] , 





wi(t,c) = Vo(t,2,0)-n + (V¢,(t,2,0)w) - n. 


We drop the higher order terms —p;¢:,(z,0,t)w and (V¢,(z,0,t)w) - 7 in these 
two equations because of the assumption of small beam displacement and obtain first 
order approximations for —p;¢; and V¢, respectively. Upon approximating the space 


domain Q(t) by N = [0,a] x [0,4], the open loop model described above is given by 


Ort = c*Ad, (x,y) = UE t > 0, 


Vo-n = 0, (2,y) €T,t>9, 
dy(t,z,0) = —w,(t,z), O0O<2<a,t>0, 
pow + 02(Elwer + cpl wert) = —psdi(t,2,0),0<a<a,t>Q0, (11.8) 
w(t,0) = w,(t,0) = w(t,a) = w,(t,a) = 0, t>0, 


o(0, x, y) do(z,y) , w(0,x2) = wo(z), 
6:(0,2,4) = di(z,y) , we(0,2) = wi(2). 
Note: Throughout this paper, 0, denotes partial differentiation with respect to the 
variable a (e.g., 0? = 2 
System II.8 is a formal representation of the dynamics of a coupled acous- 
tic/beam structure. Computational techniques (e.g., variational methods) used to 
obtain approximate solutions to this system are based on rigorous convergence argu- 
ments, typically done in the context of variational formulations of I].8. To accomplish 


this, the state is taken to be z(t) = (¢,w) in the Hilbert space H = L () x L?(T) 
with energy product 


p ¢ 
( = | i dedw + | pe wn dy , 
w n Q Cc To 
where L’(Q) is the quotient space of L* over the constant functions. Also, we define 
the Hilbert space V = H (Q) x Hé(To) where H (Q) is the quotient space of H! 
over the constant functions and Hé(To) = {% € H*(To) : o(z) = ¥,(r) = 0 at 


10 





xz =0,a}. The energy product of V is taken as 


(2 }.(5 ) = [rVo-Védw+ | Bl ssn dy. 
V 


] 





Next we define the weak variational (i.e., sesquilinear) forms: 


fi(9,8) = | py Vo-VEdw for 4,€ € H'(O) 


( 

Bo(w,n) = . WeaNer dy for w,n € HZ(To) 

pi(w,n) = I ppwndy for w,n € L*(Io) 

p2(¢, €) = fs BE dew for 4,6 € T(Q) 

Hi (w, 7) ~[ EI WeeNee dy for w,n € H3(To) (11.9) 
( 


ua(d,6) = [pj VO-VEdw for 4,€€ H(0) 

i(w,n) = f cpl Wretee dy for w,n € H3(To) 

n(n) = | pr $t,2,0)ndy for $ € H'(Q) and n € H3 (Lo) 
r2(w,£) = [ _ pr &(t,2,0)wdy for w € H3(Po) and € € A'(M) 





and express system II.8 in weak form as, 


pi(wee, 7) + Ki(w, 7) + p(w, 7) 
p2( Pr, €) + p2(¢, €) 


—T1(¢:,7) and (11.10) 
To( we, €) . (II.11) 


Our next task is to write equations II.10 and II.11 as a single second order 


differential equation. Let ® = (¢,w) and W = (£,n), such that 6, € V, and define 





sesquilinear forms: 


o1(®, WV) 


potin = | psVb-VEdw + J Elweanee dy, 
o2(®, V) = K+ —7T2 = | {ep]WeeNzz za ps( P(t, x, 0)n — E(t, c,0)w} dy , 


where 01,02 € V x V —+C (space of complex numbers). The formulations o; and 


2 satisfy coercivity and continuity (i-e., boundedness) conditions 


Roi(O,%) > ally, 


I 


A 


oi(®, Vy] < cal[Oliv|[Wllv, 
Ro(O,®) > cs(Wee; Wee) 12(Fo) = Callwllzecrg) » 


loo(®, Y) < callOllv Vly, 


where # denotes the “real part of”. The second order open loop problem is given by 
(u(t), V)vev + oa(z(t), ¥) + on(z(t), ¥) = 0, (11.12) 


where V* is the dual of space V and (-,-)y*,y denotes the usual duality pairing. 
Since o, and a» satisfy the continuity and coercivity conditions shown above, 

the Lax-Milgram Theorem guarantees the existence of uniquely determined bounded 

linear operators A,, A such that the weak and strong formulations of the coupled 


system are equivalent [Ref. 8]. That is 
(A1®, V)yxy = 01(®,V) and (42%, V)y.y = 02(0, WV). 


Thus system IJ.12 gives rise to the equivalent system defined in terms of functionals 
A, and A> 

Zu(t) + Agz(t) + Aiz(t) = 0. (11.13) 

To facilitate our numerical work, we must express system II.8 as a first order 


system. Our goal is to write Model I as 


Pt 0 0 I 0 
E | (11.14) 
Pt cA 0 0 0 od; 
Wee 0 —-Sfaf -l —2tat | | w, 
Nee pine \ cS _etrnemninsssnenssscnseensnsmnsease 
Uz A U 


The symbol II appearing in matrix A above represents ~ $,(t, 2, 0) whenever the 
product Au is calculated. 

We now develop the necessary notation and sesquilinear forms in order to write 
system II.8 as a first order system. The following approach replicates that found in 


[Ref. 9]. 
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We begin by defining the product spaces V = V x V and H = V x H whose 


norms are given by 
H(®, Y)ly = Olly + WI, and (OY) ij = [Ol + WIZ,, respectively. 


For x = (®, VW) and O = (T, A), the sesquilinear form o : Vx VY —> C is then defined 
by 

0(0, x) = o((T,A), (8, ¥)) = —(A, ®)y + 01(T,W) + 09(A,¥). 
Since the duality product (-,-)y+,v is the unique extension by continuity of the scalar 
product (-,-)y from H x V to V* x V, it follows that for appropriate restrictions on 


© we can write 


o(0,x) = o((T,A),(®,¥)) = (A, ®)y + (ArT, Vyey + (AoA, Wey 
= —(A,®)y + (AiT + AcA, V)x 
= ((—A, AT + A2A), (8, Y))x 
= (-A®O,x)x- 


The operator A:# —> H is given by 


0 I 
AL 73 


A= (1.15) 








where the domain of A= {0 = (Y,A)€ H: A EV, A\Y+A,A € H}, A, and Ay 
are the operators defined by o, and a2, respectively, and the above calculations hold 
for © € domain A. 

By letting Z(t) = (z(t), 2(¢)) and taking y € YV, the first order system is 


written in weak form as 


(Z:(t), x)vev = —(Z(t), x), 


which is formally equivalent to the strong formulation 
Z(t) = AZ(t) (11.16) 


13 





in H, where A is given in II.15. Equation IIJ.16 concisely represents the matrix 


representation of Model I given by equation II.14 where 


0 
0 0 0 I 

A= 
c7A 0 0 0 








The linear operator A is dissipative in the sense (Ay, v)x < 0 for y € domain A. To 
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see this let y = (¢, w, d:, wz), then 


Pt p 
Wr W 
(Ax; X)H = ; 
c* Ag Pt 
— 22 ow — Id, = spt Aw, Wy 


H=V xH 


= (p¢V br, VO) 12(9) + (El weet, Wee) 12(P9) + (ps AG, $2) 12/9) 


= (EI wee, Waxt) L2(Pp) —(pF Dy (t, ZL, 0), W+)12(Lo) — (cpl Wert; Wet) L2(Fo) 
Se re” LE. a earara 


by integration by parts by integration by parts 


= (pt V bs, VG)r20a)+ (pp AG, b:)2200) —- — (ps Gilt, 2,0), we) r2(75) 
\neeeernetmm, eoeenanas” 
apply Green’s formula 


(11.17) 
_ (cpl weet, Weert) 12 (To) 


= (p5Vb:, VO)12(a) + (Ps Ond, be)2(00) — (or VO, Vr) 2(0) 


_ (Pf Ot d, L, 0), We) 2(Po) = (cpl weet, Wet) 12 (To) 


(PF Ond, Pt) 12(Fo) = (Pf Pt (t, Ly 0), Wt) 12 (To) oe (cpl Weet; West) L2(To) 
qT, -oomweseorennenae” 
apply dnd=—w; on To 


(Pf Wey Pt)12 (Lo) — (Py Pelt, & 0), We) z2(r,) — (CDI Weet, Weat)12(Ty) 


= —(ceplWezt, Weat) 12 (Lo) : 


Because —cp/ || Wezel 72 (05) < 0, we see that A is in fact dissipative. 
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B. MATHEMATICAL MODEL II 





Force 


Figure 3. 2-D Acoustic Chamber for Model I] 


Our second model is similar to the first in that we consider the open loop case 
in which the wave equation is coupled with the beam equation on one side of a 2-D 
rectangular, air filled cavity as shown in Figure 3. Let 2 = [0,a] x [0,2], OX =ToUT 
with T=, Ul, UT. Let the velocity potential in 0 be given by ¢ = (t,x, y) and 
the transverse beam displacement be denoted w = w(t, x) as in Model I. Once again 
we take our perturbable boundary to be a fixed-end (i.e., clamped) Euler Bernoulli 
beam with Kelvin-Voigt damping. However, instead of assuming hard wall boundary 
conditions along I‘ as we did in the first model, we now consider dissipative boundary 


conditions along I’. Specifically, the boundary conditions are taken to be 
Ong + ad, = 0 for (z,y) ET and 0,¢= w,(t,z) for (7,0) ET, t>0, (1-18) 


where a is a constant of proportionality. As in Model I, we take f(t,z) = 0 fort > 0. 


All other assumptions remain as stated in Model I. Hence the coupled system for 
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Model I] is given by 


Dey = c27Ad, (x,y) Ee Q,t>0, 


Ond = adi, (2, y)ET,t>0, 

One = w(t,z), (z,0) ET), t>0, 
pow + O2(Elwert+ cpl wert) = —psdi(t,z,0),0<2<a,t>O, (II.19) 
w(t,0) = w,z(t,0) = w(t,a) = w(t,a) = 0, t>0, 


$(0, z, y) = go(z, y) ) w(0,2) = wo(z), 
$1(0, z, y) a 1 (2, y) } w+(0,2) = wy(z). 


The Hilbert spaces H = T’(Q) x L?(To) and V = H (9) x Hé(To) remain 
as defined in the previous section. 


In terms of the sesquilinear forms given by II.9 and 
(4,6) = [apy d€dy for 4,€€T°(), 
the weak variational form of II.19 is given by 


pi(we,7) + Ki(wi,7) + pi(w,n) = —T1(¢:,7) and (11.20) 
pal dit, €) + Ar (Gs, €) + p2(, €) T2( we, €) - (11.21) 


Here we pause to explain how the boundary condition 0,¢+ ad; = 0 gives rise to the 
sesquilinear form Ai(¢, €). Consider the wave equation ¢; = c’A¢d. By multiplying 


through by an appropriately chosen trial function € and integrating over 2 we obtain 


[  buédw = [ py AdE dw. 


Looking just at the right hand side of this equation, we use Green’s formula to obtain 


[er Ade do =f ps (Ind) dio — | py Vb VEdw where 02 =TUT». 


The boundary term represents 


a Pf (On@)§ dw = [ pr wig dy + | ps (—a) bf dy = To(wi, €) — Ar (dr, €). 
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The only difference between the weak formulations for Models I and II (equa- 
tions 11.10, [1.11 and I1.20, [1.21 respectively) is the appearance of A1(¢, &) in equa- 
tion [].21 which, as we have seen, arises because of the damping along I’. Just as the 
sesquilinear forms were shown to satisfy certain coercivity and continuity conditions 
in the previous section, so too does A1(¢, €). Once again, the Lax-Milgram Theorem 
assures us that these weak formulations give rise to uniquely determined bounded 
linear operators. Hence Model II can be expressed as a single first order equation 
analogous to equation IJ.16 in the previous section. 

Letting ® = (¢,w) and V = (€,7), such that 0, U C V, and define sesquilinear 


forms: 


o1(®, UW) 


[ o7Vb- VE dw + | Te oe 
Q To 
02(0,U) = [ {cpl WeaNez + ps((t,x,0)n — E(t, z,0)w)}dy + [aes pf ay, 


where 01,02 € V x V —+C (space of complex numbers). 


Once again, the formulations 0, and o2 satisfy coercivity and continuity con- 


ditions 
Roi($,6) > alloy, 
lor(®,Y)|| < col[@llvVllv, 
Ro(0,) 2 €3(Wee, Wee) L2(To) “Ac c3l|b|lZ2qr) = callwllzacr,) a call Pllzacry , 
lo2(®,Y) || < cal[@llv [Vly 


Denoting our state variables by z(t) = (¢,w) and making use of o, and a2 as 


defined above, we can express the second order open loop problem concisely as 
(Z(t), U)vev + o2(2(¢), V) + o1(2(t), V) = 0, (IT.22) 


where V* is the dual of space V. 
Associated with o, and a2 are functionals A,, Az such that the weak and strong 


formulations of the coupled system are equivalent. That is 


(Ai ®, W)yey = o1(®, WV) and (Ao®, W)yxy = o2(®, W) : 
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We conclude that the system described by equation II.22 gives rise to an equivalent 


system defined in terms of operators A, and A>. That is 
Z(t) + Aoz;(t) + A,z(t) = {). (11.23) 


The first order system is obtained just as it is for Model I. We let Z(t) = 
(z(t), z:(t)) and take y € V to obtain 


(Z(t), x)vev = —o(Z(t), x). 
Formally, this system is equivalent to the strong formulation 
Z,(t) = AZ(t) (11.24) 


in H. In matrix form, the strong formulation Model II is 


om 0 o 7 0 d 
We 0 0 0 I w 

= ; (11.25) 
Pet remy aN 0 0 0 Pt 
Wt 0 toe —I] —£n! Wt 
Ut A U 








—A 0 0 0 
A, = and A», = 
0 #98 298 
Po # po 


Since energy is lost along all four sides of the cavity in Model II, it is actually 
more dissipative than Model I, where energy is lost only along one side. Following 


the procedure shown in II.17 one finds that 


(Ax, X)H = —cp1||Weeel|z2(r9)—a (PF Pe, Di) L2(r) = —ep1||wellzz2@,)—aP sll bellzecry <0, 


where a is the proportionality constant which appears in boundary condition II.18. 
Now that we have variational formulations for Models I and II, we turn our 
attention to a numerical scheme often used to obtain approximate solutions to systems 


of differential equations. 
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IIT. GALERKIN APPROXIMATION 
METHOD 


In order to examine the exponential stability of approximate solutions to Mod- 
els J and H, we employ the Galerkin method to transform the systems of partial dif- 
ferential equations into systems of ordinary differential equations. A review of this 
popular computational technique follows. 

The Galerkin method is one of the variational methods (e.g., Rayleigh-Ritz, 
Least Squares, Galerkin, Collocation, etc.) which all seek good approximate solutions 
to many types of boundary-value and initial-boundary value problems of the form 
Au = f from a finite dimensional subspace. 

Specifically, when given the differential equation Au = f, where A maps the 
normed linear space X with norm ||- ||x into the normed linear space Y with norm 
|| - ly and a finite dimensional subspace Xv = span{¢y, do, ..., dv} of X, variational 


methods seek functions 
UN = C1, + CoG. +... + enOGN 
belonging to Xy which minimize 
[Aun — f lly + lluw — ullx. 


The Galerkin method is a widely used technique, subsuming both the finite el- 
ement method and the method of least squares. It is worth noting that the Galerkin 
and Rayleigh-Ritz methods coincide whenever the differential operator A is linear, 
positive definite and self-adjoint. However, there are many important differential Op- 
erators which are either nonlinear, not self-adjoint or non-symmetric for which the 
Galerkin method is applicable while the Rayleigh-Ritz method fails. For example, 
the Galerkin scheme is applicable to many parabolic and hyperbolic differential equa- 
tions whereas the Rayleigh-Ritz method may not be since the associated variational 


problem may not have a solution [Ref. 10]. 
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The Galerkin method yields a finite system of equations, from a differential 
equation, by discretizing the solution space rather than by discretizing the domain 
and operators, as is done in other popular methods. Approximate solutions are found 
using a variant of the method of weighted residuals (MWR), where the weighting func- 
tions are chosen to yield solutions which are a finite combination of known functions. 
The task at hand is to determine approximate solutions to the linear (or nonlinear) 


differential operator equation 


Au =f 


from a finite dimensional subspace Xy of some inner product space X in which 
the operator A is defined. Let (-,-) denote an inner product on X, let Xn = 
span{o, b,..., dN}, and let Yy = span{d, b’,..., bN} where Xv, Yn are N-dimen- 
sional subspaces of X. The MWR seeks an approximate solution uy to Au = f such 


that un € Xn satisfies the system of equations 
(Aun — f,o;’) = 0, forj = 1,2,...,.N. (III.1) 


Equation IJI.1 is the inner product of the residual (Auy — f) and the weighting 
function py integrated over the appropriate region. Now take uy to be a linear 


combination of the basis functions ¢’ such that 
un = Crh; + ed)’ +... + endn. 


Upon substituting this expression for uy into equation III.1, we find that ua must 
satisfy the linear system 


N | 
S> (AdMob ei =e 2 wv), for 7 = ley ere ae 
t= 1 


In the Galerkin scheme, the weighting functions pn are taken to be the basis 


functions of the approximate solution. That is 


oe = ,” for j == ee mere ke 


Ze 

















The approximate solution is obtained by solving the resulting system for the coefficient 
functions c;(t). 
An example will clarify the general procedure. Consider the following initial- 


boundary value problem 


ut — kuze = 0, overl =[0,1], t>0, (III.2) 
u(t,0) = 0, u(t,1) = 0, (III.3) 
Oe) =a: (IIT.4) 


describing the heat conduction through a thin one-dimensional rod with no sources 
of thermal energy. The temperature is held constant at one end of the rod (i.e., at 
x = 0); the other end is insulated. Here we take the thermal diffusivity of the rod, k, 
to equal one. 


In operator notation, equation III.2 can be written as 
Ut = Lu 


where the operator L = kd? and, for simplicity, we take k = 1. 

Our first task is to select appropriate basis functions which (i) possess the 
smoothness requirements of the second order problem (i.e., ¢% € C*(I)) , (ii) are 
linearly independent on I’, and (iii) satisfy the boundary conditions III.3. Here we 


choose 





For this example we take N = 2 obtaining 


N x” N : 
Py ee a and @; = >-Z, 


and because we will need the first and second derivatives of these basis functions 


later, we calculate them now: 


0.) =xz-1, Oo =1, 
8,65 =2z?—1, 626% = 22. 
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Our approximate solution takes the form 


N=2 


u(t.) ~un(t,2) = D7 ei(t)oh (2) = ea(t)oy (x) + a(t) dz (2), 


2=1 
where we seek functions c;(t). For N = 2, we need to select two weighting functions 


W1,W2 and introduce the spacial average (i.e., inner product or weighted integral) 


= d 
(w,v) = | wody, 


such that 
(apa) = (Weluny tort S12. (IIT.5) 


Note: The spacial superscript N 1s omitted below where the meaning is clear 
and an overdot denotes differentiation in time. 
For the Galerkin scheme, the basis functions ¢; serve as the weighting func- 


tions. Hence II].5 becomes 
(o;,un) = (¢;, Lun) for? =1,2 
yielding the system ~ 
(p1,¢1¢1) + (¢1, 22) — (¢1,c105¢1) — (b1, c207¢2) = 0, 


(d2,C1¢1) + (de, C2b2) — (p2, 410241) — (do, ©2052) = 0, 


or equivalently, 
(1,1) ($1, 2) | | Cy | 7 | (p1, 0261) (¢1, 0262) a | | 0 | - (qILe) 
(2,91) (P2,¢2) | | & (2, 03¢1) (G2, 02¢2) | | 2 0 


Given the particular boundary conditions for this problem, and because the basis 


functions $; we selected satisfy these boundary conditions, we find that the inner 
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product (¢;,02¢;) is equivalent to the inner product —(0,¢;,0,¢;) (ie., L is a 


self-adjoint operator). To see this, we integrate by parts: 
1 1 
(0;, 024;) = [ $1029; dz = b:02;\o = [ 0,0; : 0,0; dz , 


where the boundary term vanishes. Thus, 


Although the basis functions ¢; were initially chosen so that they belonged to 
C?(T), this result indicates that less regular functions may suffice for this problem 
(i.e, 6; € C*(L)). This fact, which is of little use in this simple example, is emphasized 
because of its theoretical import as well as its utility in reducing the computational 
complexity of more challenging problems. 


We now calculate the inner products contained in equation IJI.6 to obtain 








‘1 : | 
360 31s | | © is | | © 0 
a aE eee ae een eee Nee ome Ne, oe 
M 2 A C 0 


Thus the Galerkin method has reduced our problem of finding approximate solutions 
to a second order partial differential equation (PDE) to that of finding approximate 
solutions to a system of first order ordinary differential equations (ODEs). A much 


simpler task indeed! Multiplying equation III.7 by M1 yields 


a] | 18.9281 8.9077 | | 0 
S ~12.9231 7.1077 és 0 |. 





or equivalently, 


é, + 18.9231c, — 5.9077c, = 0, (III.8) 
2 — 12.9231c¢, + 7.1077e. = 0. (III.9) 
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To determine c,; and cz, we first solve for cp in equation IJ].9 and then differentiate 


to obtain C2 yielding 


] 
C2 = 50077 6% + 18.9231c, ) . (III.10) 
es oes 


Next substitute these equations for cp and ¢2 into equation III.9 and simplify to obtain 
Cy + 26.0308¢c, + 58.1539c, = 0. 
This is a second order ODE wth constant coefficients which is easily solved. We find 
that the general solutions for c,,c, are 
Cy = cen 24681t 4 925.5628 


~2.4681 ae~?458!* _ 23.5628 Ge7 23-5628 


Cy 
where a, @ are constants. Substituting these equations for c,, c; into equation III.10 


we obtain 


Co = 2.7853ae77 499" — 0.78548 e779 °07™ —. (III.12) 


In order to determine a and £, we make use of the given initial condition given by 
equation III.4. For (t = 0,z = .5) and (¢ = 0,2 = 1) , we know that u(0,.5) = .5 
and u(0,1) = 1, respectively. Hence we have 
u(0,.5) ~ un(0,.5) = c1(0)¢1(.5) + co(O)g2(.5) = .5 
—1.65l5a —0.01516 = .5 


and 


| 
ped 


u(0,1) ~ un(0,1) = e1(0)¢i(1) + e2(0)¢2(1) 
~2.3569a + 0.02368 


| 


From these two equations, we determine that 


a = —0.3608 and 6 = 6.3442. 
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Therefore, our approximate solution to equation III.2 is 


un (t,x) = 3 c;(t)d;(z) , 


where 


C, = —0.3608e77-46!* + 6.3449¢e—23-5628t 
c= — 1.0049 e~?-4681¢ = 4.9897 e—28-5628¢ 
2 
L 
di = eo °* 
3 
L 
$2 = 3 1 


In summary, the Galerkin method simplified the task of finding an approximate 
solution of a PDE by recasting the problem as a system of ODEs in a finite dimensional 
space. Because this example was simply meant to illustrate the Galerkin technique, 
the solution obtained provides only a very crude approximate solution to the example 
problem. Techniques for refinement of the approximate solution typically include 
such things as increasing the number of basis functions and/or selection of different 
basis functions. The interested reader is referred to [Ref. 11, 12, 13] for additional 
information regarding solution refinement techniques. In Chapter IV, we transform 


Models I and II into systems of finite dimension using the Galerkin method. 


2¢ 
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IV. FINITE DIMENSIONAL 
APPROXIMATIONS 


In Chapter II, both strong and weak variational forms were obtained for Mod- 
els I and II. Formally, at least, we demonstrated that these two formulations were 
equivalent. In Chapter III, we saw how approximate solutions to an infinite dimen- 
sional PDE could be obtained by employing the Galerkin technique to formulate 
the problem in finite dimensional spaces. Key steps in implementing this particu- 
lar discretization technique included selection of an appropriately defined set of basis 
functions, use of these basis functions as the weighting (or trial) functions, calculation 
of the inner product(s) defined for the discretization space, and finally, integration 
over the spacial domain to transform the system of PDEs into a system of time de- 
pendent ODEs. In this Chapter we extend these ideas to coupled systems of PDEs, 
to wit, Models I and II. 

Our approach will be: (i) choose finite sets of basis functions which span the 
approximating solution spaces, (ii) express the infinite dimensional state variables 
(w(t, x), o(t,z,y)) in terms of these basis functions, and (iii) use of the weak varia- 
tional forms developed in Chapter II to obtain finite dimensional representations of 
Models I and II necessary for our numerical work. 

First, let { B? ee denote the 1-D basis functions which discretize the beam 
and let {Bf }f_, denote the 2-D basis functions which discretize the cavity. For 
the moment simply note that there are n — 1 basis functions in {B?}%j and m = 
(m,,+1)(m,y+1)—1 basis functions in {Bf'}7_, where (m,+1) and (m,+1) represent 
the number of basis functions discretizing the cavity along the z, y axes, respectively. 
In Chapter V, when we look at specific finite dimensional approximations of Models I 
and H, the reader will better appreciate why these spaces have the dimensions given. 


The basis sets {BP }iy and {Bf }f,span spaces H? and H™ where the 
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n—i 
i=] 


subscripts 6 and ¢ denote ‘beam’ and ‘cavity’. That is H? = span{B?}"—} and 
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H™ = span{B?}%_,. Denoting the combined dimension of the discretized beam and 
cavity spaces as N = m+(n-—1), the approximating beam and cavity solutions are 


given by 
nom m 
w(t, 2) — S> w(t) BP (2) and p(t, 2, y) = S~ of (t) BY (2, y). 
t= k=1 


Notice that in this form the state variables are separated into products of time and 
spatial functions. 

For application of the Galerkin scheme, elements of the basis sets {B?’}7_, and 
{B?\27} serve as weighting functions for the cavity and beam, respectively. Denoting 
the product space for the first order system as HN = H™ x H™%, restriction of the 
infinite dimensional first order systems obtained for Models I and II in Chapter II to 


the space H” x H% yields 
(Z;" (t), x)1 i —o(Z%(t), M5 


for Z(t) = (¢ (t, 2, y), w(t, 2), dN (t, 2, y), w(t, 2, y))*. Note that o is model spe- 
cific (Recall that we used o to denote the first order sesquilinear vectors for both 
Models I and IJ in Chapter II.). For y = (Bj", B?), this finite dimensional first order 


equation represents the linear system 

My" (t) = AN y(t), (IV.1) 
where 
y= (9%(8), dN(O)T and OY = (A(t), 6S (t), + AN (6), 0! (Cs wh), wh)" 


denotes the N x 1 = (m+n -— 1) approximate state vector. We use an overdot to 
denote differentiation with respect to time. 

Note: Below, and for the remainder of this paper, the index ranges are k, f= 
1,---,m and 72,7 =1,---,n—1 unless otherwise noted. 

Up to this point, everything we have discussed in this chapter applies both 


to Models I and II. Now we restrict our discussion to Model I for which equation 
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IV.1 represents the linear system shown below. The non-zero entries in M% and AN 
represent block matrices. The rows of these block matrices are generated by holding 


£ and 9 fixed while varying k and i as appropriate for the particular product being 
computed. System IV.1 is given by 





Pi (ian Be BY’) 0 0 0 dx(t) 
0 P(e BP Be) 0 0 wilt) | 
0 0 po(Sye1 BE, BP) 0 bx(t) | 
0 0 0 pi(S yer BP, BP) | | w(t) 
MN y(t) 
0 0 Bi (dopa BY, BY) 0 x(t) 
0 _ 0 0 Bo(Sia1 BP, BP) | | wil?) 
—pa (Dope BE, BY) 0 0 T2(Dn1 BP, BY) | | dx(t) 
0 wi (Di BP, BP) =n (RL, BE, BP) —mi (Oe BP, Be) | | ale) 
AN y (t) 


where is =), Bo(-, a), pil-, =) pa(-, i), pal, :), Ha(-, =) ails, yi Ta(-, )s and Ki (-, ‘) refer to 
the sesquilinear forms shown in II.9. We represent this linear system concisely as 
MN 0 oN (t) 0 MN || #2) 
0 MN | | de) 


(IV.2) 
—Ay —Ay’ 











b(t) 


with 
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The component matrices are given by 


La | = | VBR. VBP de lux | = 02? 82 B? dy 
k,é 2 459 Io 


—_— P f vie) me nae Tr Nn 
MN | = Le B; dw | Ms, |. = [ mB B, dy 
— ™m m an nm92pRn 
AN iz = [over VBP dw Ay |. = | E102; 0, By dy (IV .3) 


fax | =-f oprteopray | ax | = [ o BR (e,0)Bpay 
ae Vo kj Io 


Ay, | = i cpld; BE 0B? dy, 
449 To 


? 


where the finite dimensional products correspond to the infinite dimensional sesquilin- 
ear forms given in I]J.9. 

The finite dimensional representation of Model II is similar to that given above 
for Model I, except the matrix Ad’ contains the additional sub-matrix A4\, which arises 


because of the boundary damping along I’. For Model II, A?’ is given by 


N N 
Ag As) 


N N 
Ax» Ax» 


The component matrices for Model II include all of those specified in [V.3 as well as 
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| Ay le = | ap; BRB? dy (IV.4) 


which corresponds to the weak form 1(¢,€). The linear system for Model II is given 


by 
Bi 4 BeBe) 0 0 0 b1.(t) 
0 Bo(Sspa BP, BP) 0 0 wilt) | 
0 0 pA BE BT) 0 $x.(t) 
0 0 0 pi(diy BP, Bt) | | w,(t) 
Se irre” ee jeer” 
MN y(t) 
0 0 Pid opet By , Be) 0 $x (t) 
0 0 0 Booey BP, BP) wi(t) 
—po(So_, B™, BM) 0 A Bre). BOs ee Be) $x, (t) 
0 —w1(Dynr BP, BP) —n1 (Of. BE, BP) —«i(027! BP, BP wi (t) 
\ eaten orem” 


AN y(t) 


As mentioned previously, the rows of each block sub-matrix are generated by holding 
£ and 7 fixed while varying k and 7 as appropriate for the particular product being 
computed. 

The general form of matrices MNJ" (t) and AN Y(t) is presented below to 
help the reader better conceptualize the overall structure of the system. The block 
structure of Mj‘ and A’ is characteristic of the larger matrices M% and AN for 
both Models I and H. Also, the products represented by (-,-) in the matrices below 
correspond to those given for Model I sub-matrices Mj\, Mi}, AN}, AN, and AX, in 
1V.3. We represent these matrices as 


MN 9% (t) = 
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(BP, BY) (BP, BY) ea (Br, BD) 0 ‘es 0 $y (t) 
(BY, BY) ac (Br, BD) 0 A 0 9 (t) 
(BY, Bri) =H te (BY By =i) 0 a 0 ale 
0 as 0 0 (BE, Bt) (By, Be) | | wh 
0 0 0 (BP,B3) --- (B%4,B3) | | w(t) 
0 oe 0 0 (BT, Bho) °:- (BR_1, BR») wh _»(t) 
0 pe 0 0 (Bi, Br.) ae ( BF xx Bh1) w(t) 


0 - 0 0 (BE, BP) --- (Bay, Be) bY (2) 

0 0 0 (BP,B3) ---  (BR_,, BP) ¢4 (t) 

0 0 0 (BP Be 5): tas (BE ype) oN _i(t) 

0 vee 0 0 (BT, Bri) --- (Bhi, Br) oN (t) 
(BY, BD) a os (Br, BS) (BT, By) ++. (Bhi, BY) wy’ (t) 
(BY, Bas) a ie (Br, Br-1) (BP, Ba-2) «+> (Bri, Ba-2) w_,(t) 
(oie) nes (Beg) (Bee) Brees). ee Big Br a) w_,(t) 


Dimensions of the matrices and sub-matrices associated with Models I and II 


are: 


MN,AN : 2(m-+(n—1)) x 2(m+(n—1)) 
MN, MN,AN,AN (m+ (n—1)) x (m+ (n—1) 
MN, M,AN,AN 1 mxm 
MN MN AN. AN: (n—1)x(n—1) 
Ay: mx (n—1) Aj: (n—1)xm 
In Chapter V, we obtain specific numerical approximations of M% and A® for 
Models I and II, and examine the stability of these finite systems to gain insight into 


the stability of the infinite dimensional systems II.8 and II.19. 
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Vv. SPECIFIC APPROXIMATIONS AND 
RESULTS 


In this chapter we select specific basis functions for discretization of the beam 
and cavity, discuss development of the computer programs used to generate matrices 
of the finite systems developed in the previous chapter, and present results for specific 
approximations. 

Cubic splines and tensor products of Legendre polynomials are chosen as basis 
functions for the beam and cavity spaces, respectively (We refer to the tensor prod- 
ucts of Legendre polynomials as “tensored Legendre polynomials” throughout this 
paper.). Since these choices are by no means the only possibilities, the interested 
reader is referred to [Ref. 9] for a discussion of alternate choices as well as selection 
criteria which includes: smoothness requirements, uniform preservation of exponen- 
tial stability of approximating systems, accuracy, sparsity of system matrices, and 
ease of implementation. 

The cubic splines used as a basis for H}’ satisfy the smoothness requirements 
and are easily adapted to satisfy the clamped boundary conditions. We construct 
the set {B?}%-) by first partitioning the beam into n uniform intervals of step size 
a ae Letting B? denote the standard cubic spline corresponding to this partition, 


then the basis functions for the beam discretization are taken to be 


B? = B? — 2B" — 2B", 
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where the standard cubic splines, as defined in [Ref. 11], are given by 


(oe = Gio)"; ee [a5 ee | 

h? + 3h? (ax — 23-1) + 3h(x — 24-1)? — 3(x — 2;_1)°, if x € [2;_1, 2;] 
B(x) == ¢ AP + 3h? (ain, — 2) + 3A( zig, — 2)? — 3(2i4, —2)*, if x € [2;, 2541] 

(ting —2)°, if x © [xj4y, t542] 


0, otherwise. 


All of the basis functions B? do in fact satisfy the clamped boundary conditions 
B?(0) = 6, B?(0) = B7(a) = 0,B?(a) = 0 


for 2 = 1,2,---,n —1 as can be seen in Figure 4 below where the interval [0,1] is 
partitioned into 10 uniform subintervals. The dashed curves represent B? and B”_, 
which have compact support over three intervals. The interior splines, which have 
compact support over four intervals, are denoted by solid lines. 

Tensored Legendre polynomials are used as a basis for H”. As stated in [Ref. 
9], these polynomials produce smaller, more structured matrices than those obtained 
with linear splines or finite elements, and the natural boundary conditions along 
the cavity walls obviate modification of the basis elements to satisfy some essential 
boundary conditions. Authors of [Ref. 9] assert, however, that these benefits are not 
as critical in the 2-D case as they are in the 3-D problem where system matrices are 
considerably larger. 

The basis set of tensored Legendre polynomials is obtained by forming the 
product of transformed Legendre polynomials, denoted L?(xr) and Lily), where the 
subscripts 7,7 indicate the degree of the polynomial, from the interval [—1,1] to 
(0, a] x [0,2], respectively. The transformed polynomials are obtained by substituting 
an appropriate linear transformation for x in the recursive definition of the standard 


Legendre polynomials: 


Pouce ——[(2n Phd 02) (V.1) 
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Basis Splines (Cubic) for the Beam 


Y-axis 





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
X~axis 


Figure 4. Cubic Polynomial Splines 


with Po(z) =1, P(x) =z. 


For example, substituting + for z in definition V.1 transforms the standard 
Legendre polynomials from [—1,1] to [0,1]. Orthogonality of the Legendre polyno- 
mials is preserved under this linear transformation (see Figure 5). Recalling that the 


set of basis functions for the cavity is denoted {B7"}7_,, we define B™ as 


k=1,2,---,m 
Be (x,y) = Li(z)L;(y) for 4 i=0,1,---,me+1 , (V.2) 
j =0,1,---,m,4+1 
where we impose the condition 7+ 7 4 0 to eliminate constant functions (i.e., exclude 
Lé(x)Lo(y) = 1 for all (x, y)) ensuring the set of functions is suitable as a basis for 


the quotient space. Hence, the dimension of the cavity space is m = (mz + 1)(m, + 


1) — 1. For consistency throughout this paper and in our computational algorithms, 
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Translated Legendre Polynomials 


Y-axis 





“0 01 O02 O08 O04 O5 06 O07 =O8 O89 1 
X~axis 


Figure 5. Transformed Legendre Polynomials 


the subscript & shown in V.2 is determined by holding 7 fixed while z varies. As an 


example, the indexing scheme for m, = 2 and m, = 2 is 


Br (x,y) = LY(x)Lo(y) Be (z,y) = 15(z) L4(y) 
BY (x,y) = L5(x)Loly) Be (x,y) = Lo(x)La(y) (v3) 
B3 (x,y) = Lo(x)Li(y) Br (x,y) = LY(2)L5(y) 
Be (x,y) = Li(x)Ly(y) By (2z,y) = L5(z)L5(y)- 


Having selected basis functions for the beam and cavity spaces, we now turn 
our attention to the computation of the various component matrices associated with 
Models I and IH. All computations are done using MATLAB (MATLAB is a high- 
performance interactive software package produced by The MathWorks, Inc., for sci- 
entific and engineering numeric computation.). Programs written for our numerical 


work are found in Appendix A. 
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Matrices M\, M3‘, and AX are all computed similarly since they each involve 


integration over the cavity space. In Chapter IV we defined these matrices as 


v| = | VBR.VB? du, | N | = | *erep 
a |. I k g a M3} 7 ae k De dw , 


(V.4) 
| a = [erVB VBP de. 


The matrices are computed by a routine suggested in [Ref. 9]. For the moment, 


consider the integrand of Mj) where Bf’ = L?(x)Li(y) and BY = L3(x)L'(y): 


VB? -VBY = (0,Br,0,Br) - (0,B%, 0B) 
= (0,(L7L;), @y(L7L;)) - (Ax(L5L;), Oy(L5.L5)) 
(L0,L2, L20,L;) - (Li 02%, L28,L') 

( 


t7~“"4 


Pp? "Pp 


L,L0,L20,L%) + (LE L20,L‘6,L") . 


Because of the structure of the integrand, we make use of the tensor properties of the 
transformed Legendre basis functions to construct Mi\, MX, and AN. Orthogonality 


of the transformed polynomials reduces computational complexity since 
{ a 
[ LL, = 0 and | LiL, = 0 whenever 7#q and ip. 
First construct fundamental (m, + 1) x (m, +1) matrices M™ and K™ given by 
[Mi lip= f° LE(x)Le(2)de and [KP]p= f° O,1%(2)0.13(c) dz 
0 0 


with similar definitions for Mj” and K;” (in fact, M™ = Mj? and K™ = Kj” when 
M, = m, and {0,a] = [0,/]). By orthogonality of the transformed Legendre polyno- 
mials, matrices M™ and M?” are diagonal. The matrices MN ij and MN are formed by 
computing 


MN =M? @K™+K™@M™ and My = Mz @ My. 


The symbol @ denotes the tensor product, which we accomplish by using MATLAB’s 


kron function. The ordering in the above definition is not obvious; however, attention 
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to the indexing scheme shown in V.3 and as described in Chapter IV delineate our 
convention. MN and M3\ are obtained by removing the first row and column of MN 
and MN reflecting the deletion of the constant function from the basis set. Since we 
take p; to be constant in this paper, Aiy is trivial to compute. Matrices MN and AN 
are positive definite and symmetric—although not sparse. M2. is diagonal, positive 
definite. The program matten.m, written to generate these matrices, computes the 
transformed Legendre polynomials and their derivatives iteratively. Integration of the 
differentiated transformed Legendre polynomials is accomplished by using Gaussian 


quadrature, while the orthogonality relation for Legendre polynomials, given by 
b-—a f b-a 2 
(4) P, (t) dt = ——__———6, , 
2 [Pw (t) 2 +l 





is used to compute the integrals of the translated Legendre polynomials (6,, denotes 
the standard Kronecker delta: 6,, = 0 ifr #s, =1lifr=s). Note that wat is a 
simple transformation factor which enables one to use this exact integration formula 
for integration over an interval [a, bj. 

Since matrices MN, M8, AX), and Ajj, given by 


[ux] =f aprespay, | uy] = pspresay, 
B57 0 0 


tJ 


AX, | = | EI02BP02B? dy, | aX, | = | cpld2 BY? 62B? dy, 
oe) To 49 To 


tJ 


(V.5) 


all involve cubic spline functions or their second derivatives, they are computed sim- 
ilarly. The program myspline.m is used to generate the set of basis splines. Intrin- 
sic MATLAB functions polyder and conv are used to differentiate and compute the 
product of the cubic splines. A simple program, polyint.m , performs the neces- 
sary integration. The programs used to compute M4), M4, AN, and A2, capitalize 
on the symmetry of the spline functions (mat1222.m computes Mj}, Aj,, and A>; 
matm22.m computes M?*). Because the splines are equal to zero outside their regions 
of compact support, these four matrices become seven-banded for n > 10. All four 


are symmetric and positive definite in construction. 
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The elements of matrices A}, and A, correspond to the integrated product 
of the collapsed tensored Legendre polynomials (i.e., Bf’ = L?(x)L;(y)|y=o) and the 


cubic splines over the interval [0,a]. Recall these matrices are given by 
| As = -| p Bz (x) By" (z,0)dy and | AN = | ps BR (x, 0)B? (2) , dy. 
i,t To k,J To 


These two matrices are computed using the function a9132.m . Integration, done us- 
ing Gaussian quadrature, is accomplished using MATLAB’s intrinsic function polyval 
to evaluate polynomials at translated Gaussian knots. Note that computation of each 
element of A and Aj, actually requires three or four integrations rather than just 
one (three if splines B?~* or B”7} appear in the integrand; four for integrands involv- 
ing interior splines). This is due to the piecewise construction of the cubic splines 
over their respective regions of compact support. For example, if B? is an interior 


spline, integration over Ig is given by 


p B? (x) By (x, 0) dy (x) By"(x, 0) dy 


t 


xr—-1l 1 
[Le Br(2)Br(a,0)de+ [ prBY(e)Be(2,0) da 


z—2 


Tl 
~— 
can] 
ks 
a, 
by 
sae 


Io 


zt+l o+2 
+ f° p,B(@)BR(@,0)de+ [py BP (2) Be (2,0) de 


Although A and A, are both full matrices, their computation is simplified since 
each has a well-defined block structure. Further, since we take p; to be constant, Ady 
is precisely the negative transpose of Ad, (i.e., AS, = (—A4) ). 


The function mata41.m is used to compute 
N on m RM 
nL, = [avsB By” dy 
ao; |[ BO, BEOy av+ | BE ,DBE(c,) de 
1 2 


+ [ BE (ay) BP (au) dy} . 
3 


Fundamental matrices corresponding to integration across I},[2, and [3 are com- 


| puted and then summed to generate A‘. Because of the well-defined, block diagonal 


4] 





structure of the matrices corresponding to integration across [; and I3, these ma- 
trices are computed simultaneously. The matrix produced by integration across To, 
although not sparse, possesses symmetry along diagonals which simplifies its con- 
struction. Given the structures of these fundamental matrices, Aj is symmetric with 
a well-defined block-diagonal and symmetric off-diagonal construction. 

In light of the structures of the component matrices discussed above, we note 
that for both Models I and II: (i) the matrices Aj’ and M™ are symmetric and 
positive definite by construction, and (ii) Aj’ has symmetric and skew symmetric 
block construction (A € R”*” is said to be skew symmetric if A? = —A.). 

We are now ready to examine the stability of the approximation schemes devel- 
oped for Models I and II. For our numerical work we assume the following parameters, 


which according to [Ref. 9], are physically reasonable for a .6m by 1m cavity: 


a= .6m, the Algae, py = 1.214 





c? = 1176495, pp =1.35%2, EI = 73.96 Nm? 


sec? 





epl = .00122= 


séc 


Note: For Model II, we take the proportionality constant a = 1, where a appears in 


the boundary condition (equation I1.18) 


Onp = —ad; for(z,y)ET, t>0. 
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MODEL TI: For n =m, = m, = 6,7,---,18 the margins of stability for the open 
loop system are listed in Table J. For each n, the locations of the eigenvalues, \, 
are displayed in Figures 6, 7 and 8. Eigenvalues were obtained using MATLAB’s 
toolbox function eig(A%,M™), where matrices AY and M™ are shown in equation 
IV.1. Eigenvalues having real parts with magnitude greater than 1 are excluded from 
our plots in order to better see the distribution near the imaginary axis. 

For small n, the results obtained agree very favorably with those reported in 
[Ref. 9]. Comparison of the values shown in Table I indicates that there does not 
appear to be a uniform margin of stability between the eigenvalues and the imaginary 
axis (i.e., the data in Table J does not indicate that the maximum R(\) for the cases 
tested is converging to a limit.). This is what we expect based on the conclusions 
contained in [Ref. 5]; however, we are somewhat hesitant to report this finding since 
positive eigenvalues appear in our results for n > 17. These positive eigenvalues 
should not appear since Model I is dissipative, and therefore, all eigenvalues of the 
system should lie in the left-half complex plane (i.e., R(A) < 0). The absence of 
a clear margin of stability, as well as the appearance of positive eigenvalues, may 
represent a numerical/computational instability problem. We offer two reasons for 
our suspicions. 

e During the development of the programs used to compute the component 
matrices M, AN, we were able to delay the appearance of positive eigenvalues by 
incorporating more stable computational methods. Our early programs relied heavily 

on MATLAB’s intrinsic “poly-type” functions (e.g., polyder, conv, polyval ) and our 
simple polynomial integration program polyint.m. We modified our code so that 
integrations involving tensored Legendre polynomials—or derivatives thereof—is done 
either by using Gaussian quadrature or by well-known properties of the Legendre 
polynomials. Before these changes, we observed positive eigenvalues for n = 13. 
However, after incorporating these more stable techniques, positive eigenvalues did 


not appear until n = 17. 


43 








e We first attempted to find eigenvalues by calculating D = (M™)-1.A% and 
by using MATLAB’s eig(D). However, for values of n > 13, MATLAB returned warn- 
ings that D was near singular. We then sought to solve the generalized eigenvalue 
problem using MATLAB’s eig(/ A’, M™), seeking a more computationally reliable al- 
gorithm for the problem at hand. Note that for n < 12, eig(D) and eig(/ AX, M) 
returned very similar results. The consistency of the patterns shown in Figures 6, 7, 
and 8 (all generated using eg (-,-) ) lend confidence to our belief that eig( A”, M”) 
provides more reliable results than does eig (D) (Note that ezg(-,-) did not return 
any “near singular” warnings even when tested using very poorly conditioned ma- 
trices.). Nonetheless, computational instability may increase as n does since MAT- 
LAB’s ezg(-,-) function uses a QZ algorithm and, according to [Ref. 14], some QZ 
algorithms destroy both the symmetry and positive definiteness of the semi-definite 
pair (AY, M%). 

Finally, inspection of the eigenvalue plots appearing in Figures 6, 7 and 8 
reveals a consistency in pattern even for n = 17,18, when positive eigenvalues appear. 
Thus our computational approach does not fail catastrophically for a particular (large) 
n; rather, it degenerates as the matrix systems become increasingly singular as n 


increases. 
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MODEL IT : Results obtained for n = m, = m, = 5,6,---,20 are listed in Table II. 
Figenvalue plots are displayed in Figures 9 and 10. The MATLAB toolbox function 
eig( A’, M™) is again used. Eigenvalues having real parts with magnitude greater 
than 30 are not displayed in order to better see the distribution near the imaginary 
axs. 

The data reveals that while eigenvalues lie further away from the imaginary 
axis, as expected, given the dissipative boundary conditions assumed along I as well 
as I’, no definitive uniform margin of stability appears to exist. That is, the values 
shown in Table II are creeping towards the imaginary axis as n increases. Although 
this movement cannot be seen from Figures 9 and 10, the figures do reveal consistency 
in the eigenvalue plots as n increases. Nonetheless, this creeping phenomena may not 
be an indication that the infinite dimensional system lacks uniform exponential sta- 
bility. Rather, the problem may be related to our numerical/computational approach 
for the reasons stated above. 

For Model II, we see that (i) the maximum real part of the (non-zero) eigen- 
values lie further away from the imaginary axis in this model than they do in Model I, 
and (ii) the dimension of the approximating solution spaces can be increased without 
the appearance of positive eigenvalues, (at least up ton = m, = m, = 20—the extent 
of our testing). Thus Model II is likely the better choice for use in formulating the 


noise control problem. 
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Table I]. Model I: 









maxt®O)) 
-0.021207 

ar 
8 -0.015566 
-0.006223 

















12 -0.005533 
13 -0.005440 
14 -0.006347 
-0.006806 
16 -0.001130 
17 +0.050618 






+0.051744 





Margin between the open loop eigenvalues (\) and the imaginary 






a -]. Tease 
ji an 733959 
: -1.388610 





-1.104220 
11 -1.088381 
13 -1.076611 
15 -1.065708 | 
17 -1.055837 
19 -1.047767 











Table IJ. Model II: Margin between the open loop eigenvalues (A) and the imaginary 
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—1 -0.5 0 —1 —0.5 0 
max{Re(ev)}=—0.021207 max{Re(ev)}=--0.019998 


x 10° n=mx=my=8 x 10° n=mx=my=9 





4 
~1 ~0.5 0 —1 —0.5 0 
max{Re(ev)}=-0.015566 max{Re(ev)}=—0.010471 
x 10° n=mMmx=my=10 x 10° n=mx=my=11 





—0.5 —0.5 
max{Re(ev)}=—0.006223 max{Re(ev)}=—0.004555 
x 10* n=mMmx=my=12 x 10° n=mMmx=my=13 





—1 —0.5 ) —1 —0.5 0 
max{Re(ev)}=—0.005533 max{Re(ev)}=—0.005440 


Figure 6. Mod I: Eigenvalues for n = m, = m, = 6,7,8, 9, 10,11, 12, 13. 
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—1 —0.5 0 —1 —0.5 0 
max{Re(ev)}=—0.006347 max{Re(ev)}=--0.006806 
x 10° n=mx=my=16 





—1 —0.5 0 
max{Re(ev)}=—0.001130 


Figure 7. Mod I: Eigenvalues for n =m, = m, = 14,15, 16. 
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—1 ~0.5 O : ~0.4 —0.2 O 
max{Re(ev)}=0.050618 max{Re(ev)}=0.050618 


x 10° n=mx=my=18 x 10° n=mx=my=18 





—1 0.5 O aaa 02 | Oo 
max{Re(ev)}=0.051744 max{Re(ev)}=0.051744 


Figure 8. Mod I: Eigenvalues for n = m, = m, = 17,18. 
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—30 —20 —10 0 
max{Re(ev)}=—1.088381 max{Re(ev)}=—1.087316 


Figure 9. Mod II: Eigenvalues for n = m, = m, = 5,6,7, 8, 9, 10, 11, 12. 
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Figure 10. Mod II: Eigenvalues for n = m, = m, = 13,14, 15, 16,17, 18,19, 20. 
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VI. CONCLUSIONS 


In this paper we investigate, by numerical approximation, the uniform expo- 
nential stability of two infinite dimensional systems developed to model the acous- 
tic/structure interaction of a fluid-filled, rectangular cavity (known to be dissipative). 
Model I assumes dissipative boundary conditions along one side of the boundary, while 
Model II assumes dissipation boundary conditions along all four sides of the cavity. 
We formally obtain weak variational formulations for these two models, express each 
as a finite dimensional system by discretizing the solution spaces for the acoustic pres- 
sure ¢(t, x,y) and transverse displacement of the beam w(t, x), and use the Galerkin 
technique to transform the systems of PDEs into systems of ODEs. We evaluate 
the uniform exponential stability of these systems by examining the location of their 
eigenvalues in the complex plane. Eigenvalues of these systems are determined by 


solving the generalized eigenvalue problem (AM — A’ )gj% = 0. We found that: 


e [he numerical approximations do not reflect the existence of uniform margins of 
stability for either model. The maximum real eigenvalues do not appear to be con- 
verging towards a greatest upper bound as the dimensions of the finite systems in- 
crease. Nonetheless, our numerical results clearly indicate that Model II provides a 
wider margin of stability than does Model I and, thus, is likely a better choice when 


formulating the noise control problem. 


e The choice of cubic spline and tensored Legendre polynomials—in concert with the 
use of the Galerkin method—(i) simplifies computation of the component matrices of 
M® and A%, and (ii) contributes to the overall structure of M™ and AN simplifying 


the computation of eigenvalues. 


o3 


Possibilities for future work to include: 


e Investigate alternate methods of solving the generalized eigenvalue problem rather 
than using MATLAB’s eig (AX, M). Alternate methods should attempt to capital- 
ize on the structure of the semi-definite pair (A, M). MATLAB’s eig (-,-) function 
uses the QZ algorithm and may be destroying both the symmetry and positive defi- 
niteness of the pair (A%,M™) as some QZ algorithms do. An alternate approach to 


the generalized eigenvalue problem is suggested in [Ref. 14]. 


e Assume different dissipative boundary conditions and numerically analyze the sta- 
bility of these systems using various approximation schemes. Consider models with 
medium damping and/or different coupling mechanisms between the acoustic and the 


structure components [Ref. 5]. 


e Investigate the numerical stability and preservation of exponential stability of the 
approximation schemes presented in this paper with different choices of basis functions 


to discretize the beam and cavity solution spaces, or use different schemes altogether. 


e Investigate other mathematical libraries such as NAG or IMSL, which may have 


reliable subroutines for solving the generalized eigenvalue problem presented in this 


paper. 
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APPENDIX. MATLAB FUNCTION AND 
SCRIPT FILES 


This appendix contains the programs used to compute the eigenvalues and 


produce the eigenvalue plots shown in this paper. 


1. Below are examples of the MATLAB script files which call the various function 
files shown in this appendix, as well as several intrinsic MATLAB files, for eigenvalue 


computation and plotting. 


eo 2 KK 2K 2 2K KE 2K 2K OK ok Model I, Eigenvalue Computation eK 2 2 KK 2 ok ok Bo KK 
n=6; mx=n; my=n; a=0; b=.6; el=1; rhof=1.21; rhob=1.35; c=sqrt(117649) ; 
EI=73.96; cdI=.001; [(M11,M21,A11]=matten(mx,my,b,el,rhof,c); 
[M12,A12,A22]=mat1222(a,b,n,EI,cdI); [M22]=matm22(a,b,n,rhob); 

[A31 ,A32]=a3132(a,b,n,mx,my,rhof); m=(mx+1)*(my+1)-1; nmi=n-1; t=m+nm1; 
T=zeros(t,t); M1=[Mi1 zeros(m,nm1); zeros(nmi1,m) M12]; 

M2=[(M21 zeros(m,nmi); zeros(nmi,m) M22]; M=[Mi T; T M2]; 

Ai=(-A11 zeros(m,nm1); zeros(nmi,m) -A12]; 

_ A2=[zeros(m,m) -A31; -A32 -A22]; A=[T Mi; Ai AQ]; 
ev6=eig(full(A),full(M));  e6=max(real(ev6)); subplot(2,2,1); tt=4*10°4; 
axis([-1 0 -tt tt]); hold; w=[0 0]; g=[-tt tt]; plot(w,q); 

ww=[-1 0]; qq=[tt tt]; plot(ww,qq); plot(ev6,’o’); title(’n=mx=my=6’ ) 


FA IR aR aR a ak Model II, Eigenvalue Computation So OK IR a IE fk 
n=6; mx=n; my=n; a-O; b=.6; el=1; rhof=1.21; rhob=1.35; c=sqrt (117649) ; 
EI=73.96; cdI=.001; [M11,M21,A11]=matten(mx,my,b,el,rhof,c); 
[M12,A12,A22]=mat1222(a,b,n,EI,cdI); [M22]=matm22(a,b,n,rhob); 
[A31,A32]=a3132(a,b,n,mx,my,rhof); A41]=a41(b,el,mx,my,rhof,c); 
m=(mxt+1)*(myt+i)-1;  nmi=n-1; t=mt+nmi; =zeros(t,t); 

M1i=[M1i1 zeros(m,nmi); zeros(nmi,m) M12]; 

M2=[M21 zeros(m,nmi); zeros(nmi,m) M22]; M=[M1T; T M2]; 

Ai=[-A11 zeros(m,nm1); zeros(nmi,m) -A12]; A2=[A41 -A31; -A32 -A22]: 
A=(T M1; Al A2]; ev6=eig(full(A),full(M)); e6=max(real(ev6)); 

subplot (2,2,2); tt=2*10°4; axis([-30 0 -tt tt]); hold; w=[0 0]; 

g=[-tt tt]; plot (w,q); ww=[-30 0]; qq=[tt tt];  plot(ww,qq); 
plot(ev6,’o’); title(’n=mx=my=6’) 
GAARA ORIG OK OGIO a I RIK oaks a ak oak kak ak ook 


15) 





2. Function file matten.m computes MN, MN and AX. 


OGG a CGR a ak iia ai ak 33k ak ak ak ak kak ack ak kak kak 
function [Mi1,M21,A11]=matten(mx,my,b,el,rhof,c) 


% (M11 M21 Ai1i] = matten(mx,my,b,el,rhof,c) 


4 This function produces matrices Mi1, M21, and Aii--all are 

hh (mx+1)*(nyt+1)-1 by (mx+1)*(my+1)-1 

2 Input: mx = highest degree of Legendre basis poly for x-axis 
4% my = highest degree of Legendre basis poly for y-axis 

4, b= right end point along x-axis (i.e., [0,b]) 

4 el = right end point along y-axis (i.e., [0,el]) 

7%, vhof = uniform density of fluid 

4%  c = speed of acoustic wave in fluid 

4 Written by Major J. M. Shehan, last update 21 May 95. 


7% Begin matten.m 


44% Compute Mii 
if mx == my & b == el 
mxi=mx+1; x=ones(1,mx1); vx=1:2:2*mxi; intPPx=b*(x./vx); Ma=intPPx; 
4, Determine Gaussian weights (w(i)) & evaluation points (x(i)). 
x=1:i:mx-1; x=x./sqrt((2*x+1).*(2*x-1)); j=diag(x,1)+diag(x,-1) ; 
lu xJ=eig(j); x=diag(x); [x il=sort(x); u=u(:,i); w=u(i,:).72; 
w=w’?.*2; dx=b/2; 
for i=1:mx; s=x(i); p(i,1)=1; p(i,2)=s; 
dpn(i,1)=0; dpn(i,2)=1; 
for j=2:mx 
pdi,jti)=((2*j-1) *s*p (i, j)-G-1) *pli,j-1))/j; 
dpn (i, j+1)=((2*j-1)*s*dpn (i, j)-(j-1) *dpn (i, j-1)+(2*j-1)*p(i, j))/j; 
end 
end; DPxval=dpn’ ; 
for i=1:mxl 


for j=i:mxi 
Ka(i,j)=(2/b)*sum(w’ .*DPxval (i,:) .*DPxval(j,:)); 
end 
end; 


hhhhhh Mel=Ma and Kel=Ka when b=el and mx=my. 
Mad=diag(Ma); sMad=sparse(Mad); sKa=sparse(Ka) ; 
M11T=kron(sMad,sKa)+kron(sKa,sMad); [row,col]=size(M11T) ; 
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M11=M11T(2:row,2:col); A1t=rhof*sparse(M11); %M11=full(M11); 
hhh Compute M21 

M21T=kron(Ma,Ma); t=length(M21T); 

M21=sparse(diag( (rhof/c~2)*M21T(2:t))); 


else 


Lh Compute integrals of Legendre poly’s. 
mxi=mx+1; x=ones(1,mx1); vx=1:2:2*mx1; intPPx=b*(x./vx); Ma=intPPx; 
Mad=diag (Ma) ; 
myl=my+1; y=ones(i,my1); vy=1:2:2*my1; intPPy=el*(y./vy); Mel=intPPy; 
Meld=diag (Mel) ; 
4 Compute deriv’s & eval ’product’ integrals of translated Legendre’s. 
“4 For x-axis: Determine Gauss weights (w(i)) & eval points (x(i)). 
x=1:i:mx-1;  x=x./sqrt((2*xt+1).*(2*x-1)); j=diag(x,1)+diag(x,-1); 
[u xJ=eig(j); x=diag(x); [x iJ=sort(x); u=u(:,i); w=u(1,:).72: 
w=w.*2; dx=b/2; 
for i=i:mx; s=x(i); p(i,1)=1; p(i,2)=s; 
dpn(i,1)=0; dpn(i,2)=1; 
for j=2:mx 
pt sj+1)=((2*j-1)*s*pG54)-G-1) *pCisj-1))/ i: 
dpn (i, j+1)=((2*j-1)#s*dpn (i, j)-(j-1) *dpn (i, j-1)+(2*j-1) *p (i, j))/}; 
end 
end; DPxval=dpn’ ; 
for 1=i:mxi 
for j=1i:mxil 
Ka(i,j)=(2/b) *sum(w.*DPxval (i, :) .*DPxval(j,:)); 
end 
end 


% For y-axis: Determine Gauss weights (w(i)) & eval points (y(i)). 


y=t:i:my-1; y=y./sqrt((2*y+1) .*(2*y-1)); j=diag(y,1)+diag(y,-1); 
[u yJ=eig(j); y=diag(y); [Ly ilJ=sort(y); u=u(:,i); 
wy=u(1,:).°2; wy=wy.*2; dy=el/2; 
for i=1i:my; sy=y(i); py(i,1)=1; py(i,2)=sy; 
dpny(i,1)=0; dpny(i,2)=1; 
for j=2:my 
py (i, j+1)=((2+j-1) *sy*py (i, j)-(j-1) *py (i, j-1))/j; 
dpny (i, j+1)=((2*j-1) *sy*dpny (i, j)-(j-1)*dpny (i, j-1)+(2*j-1)*py(i,j))/j; 
end 
end; DPyval=dpny’ ; 
for i=1:myl 
for j=1:nyl 
Kel (i, j)=(2/el)*sum(wy .*DPyval (i, :).*DPyval(j,:)); 
end 


of 


end 
44% Compute Mi1 & Ail 
sMad=sparse(Mad); sKa=sparse(Ka); sMeld=sparse(Meld) ; 
sKel=sparse (Kel) ; 
M1iT=kron(sMeld,sKa)+kron(sKel , sMad) ; [row,col]=size(M11T); 
M11=M11T(2:row,2:col); Aii=rhof*Mi1; % Mi1i=full(M11); 
hhh Compute M21 
M21T=kron(Mel ,Ma); t=length(M2iT) ; 
M21=sparse(diag( (rhof/c*2)*M21T(2:t))); 
end f End ’if‘ statement. 


end % End matten.m 
FIA A OO GGG GGG kok kak kaka ak kok ak ak ak a ak ak ask ak ok ak ak ok 
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3. Function file mat1222.m computes Mj}, AX, and AX. 


OC IO RGR CERIO I IIR Sk I ak ak a a ok 2k a ak a ak ake ak ak a a ake ak a ak 
function [M12,A12,A22] = mati222(a,b,n,EI,cdI) 


% (M12 A12 A22] = mat1222(a,b,n,EI,cdI) 


4 Returns M12, Ai2, and A22 matrices. 
% Input: [a,b] = domain; 


yf nh = no. of symmetric partitions of interval [a,b]; 
yA EI = stiffness coefficient; 
h cdI = damping coefficient. 


4 Note: n >= 4 required. 
4 Extrinsic functions called: mvyspline.m 
ysp 


4 Written by Major J. M. Shehan, updated 11 April 95. 


% Begin mati222.m 

4 Compute step size ’h‘ and generate ’x‘ vector. 
h=(b-a)/n; x=a:h:b; 

% Compute the cubic spline basis set for the beam. 
[B,B1,Bnmi]J=myspline(a,b,n); bi=Bi(1,:); b2=B1(2,:); b3=Bi(3,:); 

bnmii=Bnmi(1,:); bnm12=Bnmi(2,:); bnm13=Bnm1(3,:); 

4 Compute 2d derivative of cubic splines. 
ddbi=polyder(polyder(bi)); ddb2=polyder(polyder(b2)); 
ddb3=polyder(polyder(b3)); ddbnm1i=polyder(polyder(bnmii)); 
ddbnmi2=polyder(polyder(bnm12)); ddbnmi3=polyder(polyder (bnmi3)) ; 

[uu vv]=size(B); : 
for i=13:uu-12 % i=13 is index of B(2(1)) 
D2B(i-12,:)=polyder(polyder(B(i,:))); 
end 
hhh Compute M12 matrix 
ifn<4 
error(’n >= 4 required’) 

elseif n==4 
miii=polyint (conv (ddbi,ddbi) ,x(1) ,x(2)); 
m1i2=polyint (conv(ddb2,ddb2) ,x(2) ,x(3)); 
mi13=polyint (conv(ddb3,ddb3) ,x(3) ,x(4)); m11=m111+m112+m113: 


mi21=polyint (conv(ddbi1 ,D2B(1,:)),x(41),x(2)); 
m122=polyint (conv (ddb2,D2B(2,:)),x(2),x(3)); 
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mi23=polyint (conv (ddb3 ,D2B(3,:)),x(3),x(4)); m12=m121+m122+mi23; 


m131=polyint (conv (ddb2,ddbnmi1) ,x(2) ,x(3)); 
mi32=polyint (conv (ddb3 ,ddbnm12) ,x(3) ,x(4));  mi3=m131+m132; 


m221=polyint (conv(D2B(1,:),D2B(1,:)),x(1) ,x(2)); 
m222=polyint (conv(D2B(2,:),D2B(2,:)) ,x(2),x(3)); 
m223=polyint (conv (D2B(3, :) ,D2B(3,:)) ,x(3) ,x(4)); 
m224=polyint (conv (D2B(4, :) ,D2B(4,:)) ,x(4) ,x(5)); 
m22=m221+m222+m223+m224 ; 


M12=[m11i mi2 m13; mi2 m22 mi2; mi3 mi2 mii]; 

elseif n== 
m111=polyint (conv (ddb1,ddb1) ,x(1),x(2)); 
m112=polyint (conv (ddb2,ddb2) ,x(2) ,x(3)); 
mi13=polyint (conv (ddb3,ddb3) ,x(3),x(4)); mii=mii1+m112+m113; 


mi21=polyint(conv(ddb1,D2B(1,:)),x(1),x(2)); 
m122=polyint (conv (ddb2,D2B(2,:)),x(2),x(3)); 
m123=polyint (conv (ddb3 ,D2B(3,:)),x(3),x(4)); | m12=m121+m122+m123; 


mi31=polyint (conv (ddb2,D2B(5,:)),x(2),x(3)); 
m132=polyint (conv (ddb3,D2B(6,:)),x(3),x(4)); m13=m131+m132; 


m14=polyint (conv (ddb3,ddbnmi1) ,x(3) ,x(4)); 


m221=polyint (conv (D2B(1,:),D2B(1,:)),x(1) ,x(2)); 
m222=polyint (conv(D2B(2,:),D2B(2,:)) ,x(2) ,x(3)); 
m223=polyint (conv (D2B(3,:),D2B(3,:)),x(3) ,x(4)); 
m224=polyint (conv (D2B(4,:) ,D2B(4,:)) ,x(4),x(5)); 
m22=m221+m222+m2234+m224 ; 


m231=polyint (conv (D2B(2,:),D2B(5,:)) ,x(2) ,x(3)); 
m232=polyint (conv (D2B(3,:) ,D2B(6,:)),x(3) ,x(4)); 
m233=polyint (conv (D2B(4,:) ,D2B(7,:)) ,x(4) ,x(5)); 
m23=m23 1+m232+m233 ; 


m241=polyint (conv (D2B(3,:),ddbnm11) ,x(3) ,x(4)); 
m242=polyint (conv (D2B(4, :),ddbnmi2) ,x(4),x(5)); m24=m241+m242; 


m331=polyint (conv (D2B(5, :) ,D2B(5,:)),x(2) ,x(3)); 
m332=polyint (conv(D2B(6,:) ,D2B(6,:)),x(3),x(4)); 
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m333=polyint (conv (D2B(7, 
m334=polyint (conv (D2B(8,:) ,D2B(8,:)) ,x(5) ,x(6)); 
m33=mn331+m332+m333+m334; 


m341=polyint (conv(D2B(6, 
m342=polyint (conv (D2B(7, 
m343=polyint (conv (D2B(8, 
m34=m341+m342+m343 ; 





5) sD2ZB Us 2) 384) 3k (5) )% 


:),ddbnmii) ,x(3),x(4)); 
:) ,ddbnm12) ,x(4) ,x(5)); 
:) ,ddbnm13) ,x(5),x(6)); 


m441=polyint (conv(ddbnm1i,ddbnmii1) ,x(3) ,x(4)): 
m442=polyint (conv(ddbnm12,ddbnmi2) ,x(4) ,x(5)); 
m443=polyint (conv (ddbnm13,ddbnmi3) ,x(5) ,x(6)); 


m44=m441+m442+m443 ; 


M12=[mi1 mi2 m13 mi4;m12 m22 m23 m24; mi3 m23 m33 m34;m14 m24 m34 m44]; 


else 


mii1=polyint (conv(ddbi,ddb1) ,x(1) ,x(2)); 
mii2=polyint (conv (ddb2,ddb2) ,x(2) ,x(3)); 


mii13=polyint (conv (ddb3 ,ddb3) ,x(3) ,x(4)); 


m121=polyint (conv(ddb1,D2B(1,:)),x(1),x(2)); 
m122=polyint (conv (ddb2,D2B(2,:)),x(2),x(3)); 


m123=polyint (conv (ddb3,D2B(3,:)),x(3),x(4)); 


mi31=polyint (conv (ddb2,D2B(5,:)),x(2),x(3)); 


mi32=polyint (conv (ddb3 ,D2B(6,:)),x(3),x(4)); 


mi4=polyint (conv (ddb3,D2B(9,:)) ,x(3),x(4)); 


m221=polyint (conv (D2B(1, 
m222=polyint (conv (D2B(2, 
m223=polyint (conv (D2B(3, 
m224=polyint (conv (D2B(4, 
m22=m221+m222+m223+m224; 


m231=polyint (conv (D2B(2, 
m232=polyint (conv (D2B(3, 
m233=polyint (conv (D2B(4, 
m23=m231+m232+m233 ; 


m241=polyint (conv (D2B(3, 
m242=polyint (conv (D2B(4, 


>) ,D2B(1,:)),x(41),x(2)); 
:) ,D2B(2,:)) ,x(2) ,x(3)); 
1), D2BC332)) 3x (3) 24); 
>) ,D2B(4,:)),x(4) ,x(5)); 


[JD 2BC5.2)) ox 2) oes) 
:) ,D2B(6, :)) ,x(3) ,x(4)); 
:) ,D2B(7,:)) ,x(4) ,x(5)); 


>) 6 D2B(9 32) ) 4 (3) 7x4): 
:) ,D2B(10,:)),x(4),x(5)); 
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mii=miti+mi112+m113; 


mi2=m1i21+m122+mi23; 


mi3=m131i+m132; 


m24=m241+m242; 


if n== 
m25=0; % Since D2B(13,:) does not exitst as defined above for n=6. 
else 
m25=polyint (conv (D2B(4, :) ,D2B(13,:)),x(4),x(5)); 
end 
% Build M12 matrix: 
wi=[m11 m22*ones(1,n-3) miiJ]; w2=[mi2 m23*ones(i,n-4) mi2]; 
w3=([m13 m24*ones(1,n-5) mi3]; w4=[mi4 m25*ones(i,n-6) mi4]; 
M121=sparse (diag (w1)+diag(w2,1)+diag(w3,2)+diag(w4,3)+diag(w2,-1)); 
M122=sparse (diag (w3 ,-2)+diag(w4,-3)); 
M1i2=M121+M122; 
end 
A12=EI*M12; % Compute A12 
A22=cdI*M12; ¥ Compute A22 


end % End mati222.m 
OOO OOO BOR RIG RRR OK ak ak ak ak 3k ak ak ak ak ak 3k ak ak ak ak ak 3k 
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4. Function file matm22.m computes M3\ . 


ER ABABA A IA GR RA IG IIR aOR AIA Sa ak a ok ak ak ak aE ak ak ak ak ak ak 2k ak ak ake ak akeak 
function [M22] = matm22(a,b,n,rhob) 


4 [M22]=matm22(a,b,n,rhob) 


4 The function produces the (n-1)x(n-1) M22 matrix whose elements 
% are the integrals of rhob*(B(i)*B(j)) evaluated over the appro- 


4 priate partitions of [a,b] where i,j=1,2,...,n-1. B denotes 
4 cubic splines. 
y/ 


4 Input: a & b = boundary of beam, [a,b]; 

4% 2 = number of symmetric partitions of interval [a,b]: 
4 Yvhob = uniform mass density of beam 

% NOTE: n must be >= 4 for this function. 


4 Extrinsic functions called: myspline.m 
4 Written by Major J. M. Shehan, updated 8 April 95. 
4 Begin matm22.m 


4 Determine step size and build x vector. 
h=(b-a)/n; x=a:h:b; 
4 Compute cubic basis splines for beam; B1i=B(1) & Bnmi=B(n-1). 
(B,B1i,BnmiJ=myspline(a,b,n); b1=Bi(1,:); b2=B1(2,:); b3=B1i(3,:); 
4 Determine if ’n‘ is large enough and compute M22 matrix. 
if n<4 
error(’n >= 4 required’) 
elseif n== 
miii=polyint (conv(bi,b1) ,x(1),x(2)); 
mi12=polyint (conv (b2,b2) ,x(2) ,x(3)); 
mi13=polyint (conv (b3,b3) ,x(3),x(4));  mii=mi1i+m112+m113; 


mi21=polyint (conv (b1,B(13,:)),x(1),x(2)); 

m122=polyint (conv(b2,B(14,:)),x(2),x(3)); 

mi23=polyint (conv(b3,B(15,:)),x(3),x(4)); m12=m121+m122+m1i23: 
mi31=polyint (conv(b2,Bnmi(1,:)),x(2),x(3)); 
m132=polyint (conv (b3,Bnm1(2,:)),x(3),x(4)); m13=m131+m132; 
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m221=polyint (conv(B(13,:),B(13,:)),x(1),x(2)); 
m222=polyint (conv(B(14,:) ,B(14,:)),x(2),x(3)); 
m223=polyint(conv(B(15,:) ,B(15,:)),x(3) ,x(4)); 
m224=polyint (conv(B(16,:) ,B(16,:)),x(4) ,x(5)); 
m22=m221+m222+m223+m224 ; 


M22=rhob*[m11 m1i2 m1i3; mi2 m22 mi2; m13 mi2 mii] 
elseif n== . 
mi11=polyint(conv(bi,b1) ,x(1),x(2)); 
mi12=polyint (conv(b2,b2) ,x(2) ,x(3)); 
m113=polyint (conv(b3,b3) ,x(3),x(4)); mii=mi11i+m112+m1i13; 


m121=polyint (conv(b1,B(13,:)),x(1),x(2)); 
m122=polyint (conv(b2,B(14,:)),x(2),x(3)); 
mi23=polyint (conv(b3,B(15,:)),x(3),x(4)); m12=m121+m122+m123; 


mi31=polyint (conv(b2,B(17,:)),x(2),x(3)); 
m132=polyint (conv(b3,B(18,:)),x(3),x(4)); m13=m131+m132; 


mi4=polyint (conv(b3,Bnmi(1,:)) ,x(3) ,x(4)); 


m221=polyint (conv (B(13,:),B(13,:)),x(1) ,x(2)); 
m222=polyint (conv (B(14,:),B(14,:)),x(2),x(3)); 
m223=polyint (conv(B(1i5,:),B(15,:)),x(3) ,x(4)); 
m224=polyint (conv(B(16,:) ,B(16,:)),x(4) ,x(5)); 
m22=m221+m222+m223+m224 ; 


m231=polyint (conv (B(14,:),B(17,:)),x(2) ,x(3)); 
m232=polyint (conv(B(15,:),B(18, :)),x(3) ,x(4)); 
m233=polyint (conv(B(16,:) ,B(19,:)),x(4) ,x(5)); 
m23=m231+m232+m233 ; 


m241=polyint (conv(B(15,:),Bnm1(1,:)),x(3) ,x(4)); 
m242=polyint(conv(B(16,:) ,Bnmi(2,:)),x(4),x(5)); m24=m241+m242; 
M22=rhob*[mii mi2 mi3 m14;m12 m22 m23 m24:m1i3 m23 m22 mi2;m14 m24 m1i2 mii]; 


else 
miii=polyint(conv(b1,bi) ,x(1),x(2)); 
m112=polyint (conv(b2,b2) ,x(2),x(3)); 
mi13=polyint (conv(b3,b3) ,x(3),x(4)); mii=mi1i+mii2+m1i1i3; 
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mi21=polyint (conv(b1i,B(13,:)),x(1),x(2)); 
mi122=polyint(conv(b2,B(14,:)),x(2),x(3)); 
m123=polyint(conv(b3,B(15,:)),x(3),x(4)); m12=m121+m122+m123;: 


| 
m131=polyint (conv(b2,B(17,:)),x(2),x(3));: 
mi32=polyint (conv(b3,B(18,:)),x(3),x(4)); m13=m131+m132; 


m14=polyint (conv(b3,B(21,:)),x(3) ,x(4)); 


m221=polyint (conv(B(13,:) ,B(13,:)),x(1),x(2)); 
m222=polyint (conv(B(14,:),B(14,:)),x(2),x(3)); 
m223=polyint (conv(B(15,:),B(15,:)),x(3) ,x(4)); 
m224=polyint (conv(B(16,:),B(16,:)),x(4),x(5)); 
m22=m221+m222+m223+m224 ; 


m231=polyint (conv(B(14,:) ,B(17,:)),x(2) ,x(3)); 
m232=polyint (conv(B(15,:),B(18,:)),x(3),x(4)); 
m233=polyint (conv(B(16,:),B(19,:)),x(4),x(5)); m23=m231+m232+m233; 





m241=polyint (conv(B(15,:) ,B(21,:)),x(3) ,x(4)); 
m242=polyint (conv(B(16,:) ,B(22,:)),x(4),x(5)); m24=m241+m242; 


m25=polyint (conv (B(16,:) ,B(25,:)) ,x(4) ,x(5)); 


4 Build M22 matrix: 
wi=[mi1 m22*ones(1,n-3) mii]; w2=[m12 m23*ones(1,n-4) m12]; 
w3=[m13 m24*ones(i,n-5) m1i3]; w4=[m14 m25*ones(1,n-6) mi4]; 
M22i=diag(w1) + diag(w2,1) + diag(w3,2) + diag(w4,3); 
M222=diag(w2,-1) + diag (w3,-2) + diag(w4,-3); 
M22=M221+M222; M22=rhob*M22; 
end 
end ¥% End matm22.m 
FOE OIC AGGRO IGOR IGRI CI a ask a ak 2k ak ak ak ak ake ak ak ak 
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5. Function file a3132.m computes A} and Ad. 


ROR GIG GGG IE sk a a kak kak 
function [A31,A32]=a3132(a,b,n,mx,my ,rhof) 


% {A31 A32] = a3132(a,b,n,mx,my,rhof) 


% Currently, this function returns matrices A31 and A32. The elements 
4 of these matrices correspond to the integrated product of the 

% collapsed tensored Legendre poly’s (i.e., y=0) and the cubic poly’s 
% over La,b]. Note: The cubic poly’s satisfy clamped beam boundary 


% conditions. A32 is formed by computing -A31’. Integration is 

4 accomplished using Gaussian quadrature. 

h 

4 Input: [a,b] = interval of integration (i.e., length of beam) 

y/ n = number of symmetric partitions [a,b] is divided into 

y/ mx = highest degree of Legendre poly in basis set for beam 
i my =" ‘ " ‘ N " for cavity 

y/ rhof = density of fluid 

4 Extrensic functions called: legtrans.m to compute Legendre poly’s 
hh myspline.m to compute cubic splines 


%Z Written by Major J. M. Shehan, 13 May 95. 
4A Begin a3132.m 


7%, Compute Gaussian quadrature weights and knots for partitioned beam. 
p q g& p 
h=b/n; v=0:h:b; 


k=round((4+mx)/2); % Determine no. of knots. 
x=1:4:k-1;  x=x./sqrt((2*x+1).*(2*x-1)); j=diag(x,1)+diag(x,-1); 
lu xJ=eig(j); x=diag(x); [x i]J=sort(x); u=u(:,i); w=u(41,:).°*2; 
w=w.*2; % ?w' denotes weights; ’x‘ denotes knots. 
/, Translate knots to appropriate interval. 
dx=h/2; u=1:2:2*n; x=x’; 
for 1i=i:n 
X(i,:)=dx*(xtu(i)); 
end 
4 Compute Legendre & cubic poly’s. 
[L]=legtrans(b,mx); [B,B1,Bnmi]=myspline(a,b,n); 
% Evaluate Legendre’s at translated knots corresp to beam partitions. 
P(i:n,1:k)=ones(n,k); s=0; 
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for 1=2:mxt+1 
for j=i:n 
P(j+tnts,:)=polyval(L(i,:),X(j,:)); 
end; s=stn; 





end 
% Compute integrals involving B(1) and B(n-1) cubic splines. 
M=zeros(n-1,mxt+1); % Allocate storage space. 
qi=w.*polyval(B1(1,:),X(1,:)); q2=w.*polyval(B1(2,:),X(2,:)); 
q3=w.*polyval (B1(3,:),X(3,:)); 
ri=w.*polyval(Bnm1i(1,:),X(n-2,:)); 
r2=w.*polyval (Bnmi(2,:),X(m-1,:)); 
r3=w.*polyval (Bnm1(3,:),X(n,:)); 
s=0; 
for i=1:mx+i 
vi=sum(q1i.*P(its,:)); v2=sum(q2.*P(2t+ts,:)); 
v3=sum(qg3.*P(3+s,:)); 
ul=sum(ri.*P(n-2+s,:)); u2=sum(r2.*P(n-its,:)); 
u3=sum(r3.*P(n+s,:)); 
M(1,1)=dx*(vitv2+v3); M(n-1,i)=dx*(uitu2+u3); s=s+tn; 
end 
_f Compute interior integrals (B(i),P(j)) for i=2,...,n-2 & for j=i:mx+1. 
[uu vv]=size(B); 
s=0; r=0; 
si=w.*polyval(B(13,:),X(1,:)); s2=w.+*polyval(B(14,:),X(2,:)); 
s3=w.*polyval(B(15,:),X(3,:)); s4=w.*polyval(B(16,:),X(4,:)); 
for i=1: (uu-24)/4 
for j=i:mxtl 
ti=sum(s1.*P(i+tstr,:)); t2=sum(s2.*P(2+str,:)); 
t3=sum(s3.*P(3+str,:)); t4=sum(s4.*P(4+s+tr,:)); 
M(Citi1, j)=dx*(t1+t2+t3+t4);  s=stn; 
end; s=0; r=i+tr; 
end; M=M’; 
4 Generate A31 & A32 using M and fact that P(y=0)=(-1) or (1) for all 
Acollapsed Legendre poly’s. 
m=(mxt+1)*(my+1)-1;  A31=zeros(mti,n-1);  s=0; 
for i=1: (my+1) | 
A31(1+s:mxt+its,:)=(-1)7 (iti) *M; s=stmxt+i; 
end 
A31=-rhof*A31(2:mt1,:); A32=-A31?’; 


SS scan 








end %¥% End a3132.m 
FARO IA GR GOR IORI I aI a a ak aE ak a Ik aa ak ak ak akc a akc ak aka ake ak ake ak ak ak 
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6. Function file mata{1.m computes A . 


OG OA RAC GR aa a a kaka aK 3 a ak a ak ak ak ak akeak 6 
function [A41]=mata41(b,el,mx,my,rhof) 


4 (A41]=mata4i(b,el,mx,mny,rhof) 


4 Computes matrix A41 for Model II. 
4 Input variables: b = right bdary of [0,b]. 


hh el = upper bdary of [0,el1], 

yA mx = highest degree of Legendre poly for bean, 

h my = highest degree of Legendre poly for cavity, 
yA cc = proportionality constant of damping tern. 

4 In this program, matrix A41A corresponds to integration over 

4 0<=y<=el, x=0; A41B corresponds to integration over 0<=x<=b, y=el; 
4 and A41C corresponds to integration over 0<=y<=el, x=b. 


4 Written by Major J. M. Shehan; last update: 23 May 95. 
mxi=mx+1; myi=my+i; m=mxi*my1-1; mi=m+i; % Notation simplification. 


Ahh This algorithm can be used to computes A41C. 

4% yzones(1i,my1); v=1:2:2*myi; intPP=el*(y./v); T=ones(mxi,my1); 
4 A41C=zeros(mi,mi); q=0; g=0; 

4 for i=1:my1 

yA A41C(1+q:mxit+gq,1i+g:myl+g)=intPP(i)*T; q=mxit+tq; g=myitg; 

4% end; A41C=2+*A41C(2:m1,2:m1) ; 


4%% This algorithm computes A41A & A41C simultaneously. 
y=ones(1,my1); vy=1:2:2*myi; intPPy=el*2*(y./vy); A41AC=zeros(mi,m1) ; 
s=0; q=mx1i; g=myi; 
for i=1:mxi 
for j=1:myl 
Block (i, j)=1+(-1)* (jt+tts) ; 
end; s=st1; 
end 
for 1=2:ny1 
A41AC(itq:mxitq,itg:myl+g)=intPPy(i)*Block; g=mxi+q; g=myitg; 
end; A41AC=A41AC(2:m1i,2:m1); 
A44AC(1:mx,1:my)=intPPy(1)*Block(1:mx,2:my1) ; 
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Ahh This alogrithm computes A41B. 
x=ones(1,mx1); vx=1:2:2*mx1; intPPx=b*(x./vx); t1=intPPx(1); 
t2=intPPx(mx1); intPPx(1)=t2; intPPx(mx1)=t1; D=diag(intPPx) ; 
A41B=zeros(mi,m1); g=0; 
for i=i:myl 
A(i:mxi,i+g:mxit+g)=D; g=g+mx1;. 
end 
for k=1:my 
A41B (k* (mx1) +1: (k+1)*(mx1),:)=A; 
end; A41B(1:mx1,:)=A; A41B=A41B(1:m,1i:m); 


Ahh Compute A41 
A41=-rhof*(sparse(A41AC)+sparse(A41B)) ; 


end % End mata4i.m 
SOO GG GOO GCG RII I IRI aI AIGA IIE ak ocak ak ae ak ak ake ake a ak ak 
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7. Function file myspline.m computes the cubic splines. 


2K IK 2K 2K 26 2K 2 2K 2K 2K 2 2 2K 2 2 2 2c 2 2 2c ik 2 2c 2c 2 2 2k ok oe ice 2c 2 2 2 ko 2g 2 ok 2k 2k ok 2 2k 2 2 ke KK KK 2K KE 2 OK OE OK KE RK OK KK KE 


function [B,B1,Bnm1] = myspline(a,b,n) 

% ([B,Bi,Bnmi] = myspline(a,b,n) 

4 This function returns a family of standard cubic splines 

4 defined over [a,b], as well as the "boundary splines" which 
4 satisfy clamped boundary conditions. The interval [a,b] is 


4 divided into ’n‘ uniform partitions; ’h’ is the step size. 


% Input: [a,b] = interval along x-axis; 


y/ n = no. of equispaced partitions of [a,b]. 

% OQutput: B = each row of "matrix" B corresponds to a cubic poly 
h which is defined only over one step size (b-a)/n; 
y/ every 4 rows constitute a piecewise smooth cubic 

y/ polynomial which is non-zero only over 4 intervals 
y/ (e.g., rows 1-4 is the first cubic poly, rows 5-8 
y/ makes up the second basis function, etc.). 

h Bi = left most cubic spline satisfying clamped boundary 
y/ conditions. 

h Bnmi = right most cubic spline satisfying clamped boundary 
y/ conditions. 


4 Written by Major J. M. Shehan, updated 10 April 95. 
4 Begin myspline.m 

7, Form standard cubic splines. 

h=(b-a)/n; x=a-3*h :h:b+3*h; hh=1/h73; z=0; 


for 1 = 1:nt3 
BCitz,:)=hh*[1  -3*x(i) 3%*x(i)72 -x(i)73]; 
B(ititz,:)=hh*[-3 3*h+9*x(iti) 3*h*2-6*h*x (iti) -9¥*x(Cit1)72 
h*3-3*h* 2*x (it+1)+3*h*x (i+1) 724+3*x (iti) 73]; 
B(it2+z,:)=hh*[3 3*h-94#x(it+3) 9*x(i+3) *2-6*h*x (i+3)-3*h7*2 
h~3+3*h* 2*x (i+3) +3*h*x (i+3) ~2-3*x (i+3) 73]; 
B(it3+z,:)=hh*[-1 3%*x(it4) -3*x(it4)°2 x(it4)73]; 
Z=Zt3 ; 
end 
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“4 Form exterior splines which satisfy clamped boundary conditions. 
4 Form B1i=B(0)-2*B(1)-2*B(-1). 
bi=B(7,:)-2*B(10,:)-2*B(4,:); b2=B(8,:)-2*B(11,:): 
b3=-2*B(12,:); 
Bi=[bi; b2; b3]; 
4 Form Bnmi = B(n) - 2*B(n+1) - 2*B(n-1). 
(uu vv]=size(B) ; 
bnm11=-2*B(uu-11,:); bnm12=B(uu-7, :)-2*B(uu-10,:); 
bnmi3=B(uu-6, :)-2+*B(uu-9, :)-2*B(uu-3,:): 
Bnmi=[bnmii; bnmi2; bnmi3]; 
end % End myspline.m 
COO OEIC IC A IG Ia Ek i a ak ak 3k a ak ak ak ak feake a ak 
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8. Function file legtrans.m computes the transformed Legendre polynomials. 


CTeCret et CCT tTrtT tt rter rrr Cr rere ceotor oor efor LS eee et et Ser re re ee 
function [L] = legtrans(b,n) 


% ([L] = legtrans(b,n) 


4, This function produces a ’matrix’ L whose rows are translated 

7, Legendre polynomials of Oth through nth degree defined on [0,1]. 
4% The Oth degree polynomial corresponds to the first row of the 

4 output matrix, while the nth degree polynomial corresponds to the 
4, nti row (i.e., the last row) of the output matrix. 


4% Input arguments: n = highest degree of translated Legender poly 
hh desired; 
yA b = right end point of interval assuming [0,b]. 


% Algorithm written by Major J. M. Shehan, updated 11 April 95. 


4 Begin legtrans.m 


4 Generate tranlated Legendre poly’s of Oth & ist degree. 
L(1,:)=[zeros(1,n) 1]; L(2,:)=[zeros(i,n-1) 2/b -1i]; 


4, Generate 2d-nth deg translated Legendre poly’s recursively. 

k=0; r=nt+1; 

for 1=2:n 
d=i+1; p=[(2*(i-1)+1)*2/b (2*(i-1)+1)*(-1)]; 
L(d,:)=(1/i)*([zeros(1,n-i) conv(p,L(i,n-k:r))J - [(i-1)*L(i-1,:)]); 
k=k+1; 

end 


end % End legtrans.m 
EAE COCA IO GG AAG IR CICIGI OEIC I Ia aa AI 1 a a A Ik 1 a a AE 
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9. Function file polyint.m performs polynomial integration. 


A 2K 2 2 2 2 2 26 26 2 2 2 ok 2K 3 2K 2 6 og 2 a 2 2 ie i 2 2 2 6 2 2 2 i ok 2k 2k 2 2 i oo E22 oe 2 2 2 ok 2c ok ok 2 2 2 2 2K 2 22k ok 2 2 2 2k akc ok ok 2k 2k 


function [polyint] = polyint(p,a,b) 

4 This function integrates the polynomial ’p’ over the interval [a,b]. 

4 The polynomial ’p’ is written as a vector ’v’ with coefficients listed 
4 in descending order (e.g., 3x°2 + 5x - 8 ===> [35 8]). 

, Written by Major J. M. Shehan 10 Feb 95. 

4 Begin polyint.m 


v=[p 0]; y=v./[Llength(p):-1:1 1]; polyint = polyval(y,b) - polyval(y,a); 


end % End polyint.m 
EGOS EAHA AIGA AICI IK ak a ak ak a a 53k 3 a4 2k ak 
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